Gtx  mm 

wwBtiwaw 


* 


THE  UNIVERSITY  OF  ALBERTA 
RELEASE  FORM 

NAME  OF  AUTHOR  Arokia  Nathan 

TITLE  OF  THESIS  Numerical  Analysis  of  MOS  Magnetic  Field  Sensors 
DEGREE  FOR  WHICH  THESIS  WAS  PRESENTED  Master  of  Science 
YEAR  THIS  DEGREE  GRANTED  Fall  1984 

Permission  is  hereby  granted  to  THE  UNIVERSITY  OF  ALBERTA  LIBRARY  to 
reproduce  single  copies  of  this  thesis  and  to  lend  or  sell  such  copies  for  private,  scholarly 
or  scientific  research  purposes  only. 

The  author  reserves  other  publication  rights,  and  neither  the  thesis  nor  extensive 


extracts  from  it  may  be  printed  or  otherwise  reproduced  without  the  author’s  written 

/  // 

permission. 


CO 


THE  UNIVERSITY  OF  ALBERTA 


NUMERICAL  ANALYSIS  OF  MOS  MAGNETIC  FIELD  SENSORS 


by 


ARQK.IA  NATHAN 


A  THESIS 

SUBMITTED  TO  THE  FACULTY  OF  GRADUATE 
STUDIES  AND  RESEARCH  IN  PARTIAL  FULFILMENT 
OF  THE  REQUIREMENTS  FOR  THE  DEGREE  OF 
MASTER  OF  SCIENCE 

DEPARTMENT  OF  ELECTRICAL  ENGINEERING 


EDMONTON,  ALBERTA 


FALL,  1984 


' 


THE  UNIVERSITY  OF  ALBERTA 


FACULTY  OF  GRADUATE  STUDIES  AND  RESEARCH 


The  undersigned  certify  that  they  have  read,  and  recommend  to 
the  Faculty  of  Graduate  Studies  and  Research,  for  acceptance,  a  thesis 
entitled,  "Numerical  Analysis  of  MOS  Magnetic  Field  Sensors," 
submitted  by  Arokia  Nathan  in  partial  fulfillment  of  the  requirements 
for  the  degree  of  Master  of  Science. 


To 


My  Mother 


and  Father 


(iv) 


Abstract 


Numerical  solutions  in  two  dimensions  have  been  computed,  for  the  system  of  coupled 
nonlinear  partial  differential  equations  which  govern  drift,  diffusion,  electrostatic  potential, 
generation  and  recombination  in  a  semiconductor  device  in  the  presence  of  a  magnetic  field.  Two 
types  of  device  geometries  are  considered  ;  (a)  rectangular  silicon  slabs  under  doped  and  intrinsic 
conditions  (b)  split -drain  N -channel  MOSFET's  operating  in  the  linear  region.  The  sensitivities 
associated  with  the  various  M05FET  channel  geometries  are  studied  and  a  circuit  configuration 
incorporating  a  M05FET  magnetic  field  sensor  is  presented. 
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1.  INTRODUCTION 


A  magnetic  field  sensor  (MFS)  is  an  input  transducer  which  converts  the  magnetic 
induction  B  into  an  electronic  signal.  Many  MFS  exploit  the  Lorentz  force  F  =  q(v  x  B)  on 
electrons  in  a  metal,  a  semiconductor  or  an  insulator  in  some  manner  [1]. 

We  can  distinguish  between  direct  and  indirect  MFS  applications.  In  direct  applications, 
the  MFS  gathers  information  about  field  magnitude  and/or  direction.  These  applications  include 

H.2] 

i.  measurement  of  the  earth's  magnetic  field 

ii.  reading  heads  for  magnetic  tapes/disks 

iii.  magnetic  apparatus  control 

In  indirect  applications,  the  magnetic  field  is  used  as  an  intermediary  carrier  for  detecting 
nonmagnetic  signals.  These  applications  include  [1,2] 

i.  contactless  switching 

ii.  linear  and  angular  displacement  detection 

iii.  wattmeters 

iv.  biomagnetometry 

These  groups  each  have  their  own  sensor  requirements  since  the  variation  of  field  strengths  from 
case  to  case  is  large.  Hence,  this  broad  range  cannot  be  covered  by  only  one  type  of  sensor.  The 
sensitivity  of  the  devices  considered  in  this  thesis,  is  sufficient  for  the  indirect  applications  (i),  (ii) 
and  (iii).  In  these  applications  the  magnetic  induction  is  in  the  range  of  5-100  mtesla. 

There  are  semiconductor  magnetic  sensors  based  on  Si,  GaAs  or  InSb  technologies.  Silicon 
offers  the  advantage  of  cheap  mass  manufacturing  by  integrating  a  basic  sensor  element  with 
support  and  processing  circuitry  in  bipolar  or  CMOS  technology.  The  companies  involved  in  the 
research  and  development  of  MFS  include  IBM,  Philips,  Honeywell,  Landis  &  Gyr,  Toyota. 
Sensors  developed  outside  main -stream  IC  technologies  face  extra  cost  and  tools  as  well  as 
reliability  and  test  procedures  [1].  The  disadvantages  in  using  Si  MFS  are  its  temperature 
limitations,  e.g.  150  C  is  usually  the  quoted  limit  for  Si  device  operations.  Si  is  not  an  outstanding 
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magnetosensitive  material  as  its  mobility  values  are  inferior  to  that  of  GaAs  and  InSb.  Toshiba 
(Japan)  and  Siemens  (Europe)  produce  MFS  based  on  GaAs  technology.  GaAs  MFS  are  useful  at 
temperatures  beyond  150  C  and  GaAs  has  a  mobility  which  is  5  times  higher  than  that  of  Si. 
Research  in  magnetic  field  sensors  based  on  InSb  thin  films  or  crystal  is  strongly  pursued  in  Japan 
by  Electrotechnical  Laboratory  (ETL)  and  Hitachi.  InSb  has  an  extremely  high  mobility  (50-60 
times  higher  than  that  of  Si)  but  has  a  drawback  of  a  small  band  gap  (0.25  eV). 

The  modeling  of  semiconductor  MFS  reported  in  this  thesis  is  based  on  Si  but  sensors  of 
other  materials  like  GaAs  or  InSb  can  easily  be  modeled  by  suitably  changing  the  appropriate 

•k 

material  parameters,  such  as  p,  p,  e  etc.. 

In  the  absence  of  a  magnetic  field,  analytical  and  numerical  modeling  of  semiconductor 
devices  is  a  well  established  art  [3-5].  Many  devices  are  generally  modeled  by  the  numerical 
solution,  in  two  dimensions,  of  the  well  known  three  coupled  nonlinear  elliptic  partial  differential 
equations  comprising  the  continuity  equations  and  the  Poisson's  equation  [3-5]. 

But  numerical  modeling  of  magnetic -field -sensitive  silicon  microtransducers  has  not  kept 
up  with  the  pace  of  the  rapid  research  and  development  of  such  sensors  (see,  e.g.  [6-12]  and 
references  therein).  These  sensors  include  bipolar  and  lateral  magnetotransistors,  split-drain 
MOSFET's  and  CMOS  magnetodiodes.  The  analysis  of  these  devices  explained  in  terms  of  simple 
analytical  models  (such  as  Lorentz  deflection,  Hall  voltage,  emitter  efficiency  modulation  and 
magnetoconcentration)  has  led  to  some  dispute  [13,14].  The  analytical  modeling  of  such  devices  is 
not  straightforward  for  the  following  reasons  ; 

i.  the  basic  magnetic  field  effects  in  a  semiconductor  are  two  dimensional,  e.g.  Hall  field, 
Lorentz  deflection  and  the  magnetoconcentration.  One  or  the  other  of  these  effects  may 
prevail  in  different  pans  of  the  same  device.  An  analytical  superposition  of  the  basic  effects  is 
in  general  not  possible.  For  instance,  it  is  not  clear  which  mechanism  prevails  generally  in 
magnetotransistors  [13,14]. 

ii.  the  distributions  of  carriers  and  potential  in  a  magnetic -field-sensitive  device  are  governed  by 
the  two  interconnected  mechanisms  of  the  Lorentz  deflection  and  boundary  conditions  ;  for 
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instance,  there  is  always  a  strong  change  in  either  the  carrier  concentrations  or  the  electric 
potential  near  the  open  boundaries  of  such  a  device.  Only  if  both  magnetoconcentration  and 
space  charge  effects  are  negligible,  will  conformal  mapping  techniques  [15]  yield  a  feasible 
method  for  description  of  such  a  device. 

For  these  reasons,  numerical  modeling  of  magnetic  field  sensors  is  a  highly  desirable  tool 
towards  a  better  theoretical  understanding  of  the  electronic  mechanisms  underlying  the  operation 
of  these  devices.  The  modeling  of  magnetic -field -sensitive  semiconductor  devices  is  made 
particularly  difficult  by  the  fact  that  the  magnetic  field  breaks  some  of  the  symmetries  exploited  in 
the  modeling  of  zero  field  semiconductor  devices.  This  may  be  the  reason  why  the  numerical 
analysis  of  semiconductor  devices,  in  the  presence  of  a  magnetic  field,  was  done  only  recently 
[16,17]. 

Presented  in  this  thesis  are  numerical  solutions  in  two  dimensions  of  the  system  of  coupled 
nonlinear  partial  differential  equations  which  govern  carrier  drift,  diffusion,  electrostatic 
potential,  generation  and  recombination  in  semiconductor  devices  in  the  presence  of  a  magnetic 
field  and  under  appropriate  physical  boundary  conditions.  The  thesis  starts  with  the  discussion  of 
the  carrier  transport  equations,  in  the  presence  of  a  magnetic  field  and  the  deriving  of  appropriate 
relations  between  J  and  B.  Throughout  the  thesis  the  magnetic  field  is  considered  perpendicular 
to  the  device  plane.  A  brief  review  is  given  of  the  basic  magnetic  effects,  such  as  Hall  effect, 
carrier  deflection  and  magnetoconcentration.  Based  on  each  of  these  effects  are  devices  realised  in 
both  bipolar  and  MOS  technology  and  a  few  examples  of  such  devices  are  covered. 

Chapter  3  gives  the  analysis  of  a  homogeneous  silicon  slab,  starting  with  a  review  of  the 
basic  equations,  boundary  conditions,  and  physical  models  of  the  carrier  mobility  and 
recombination.  Then,  in  the  section  on  numerical  method,  the  discretisation  and  linearisation  of 
the  coupled  nonlinear  partial  differential  equations  as  well  as  the  boundary  conditions  at  the 
"floating"  (open)  boundaries  are  discussed  in  detail.  Here,  the  discretisation  of  the  continuity 
equations  uses  the  generalisation  of  the  Scharfetter-Gummel  scheme  to  two  dimensions  and 
nonzero  magnetic  field.  The  solution  technique  -  SLOR  (Successive  Line  Over  Relaxation),  an 
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iterative  technique  to  solve  the  linear  system  of  equations  (after  a  Newton  iterative  step)  -  is 
described  by  a  flow  chart.  The  results  of  the  slab  analysis  are  covered  in  appendices  D,  E  and  F. 
Appendices  D  and  E  show  the  influence  of  the  geometrical  aspects  of  uniformly  doped  silicon  slabs 
on  the  electric  potential  and  current  densities.  Appendix  F  shows  the  effect  of  the 
magnetoconcentration  on  the  electron  and  hole  densities  in  slabs  under  intrinsic  or  "close  to" 
intrinsic  conditions,  giving  rise  to  space  charge  effects. 

In  Chapter  4,  a  magnetic  field  sensor  based  on  the  split -drain  configuration  of  the 
MOSFET  is  analysed  for  the  linear  region  of  operation.  The  analysis  first  presents  the  basic 
assumptions  which  reduce  the  problem  (which  is  originally  three  dimensional)  to  two  dimensions 
by  integrating  the  current  density  over  the  depth  of  the  channel.  This  is  followed  by  a  description 
of  the  basic  equations  and  boundary  conditions.  The  mobility  model  employed  makes  use  of  a 
decomposition  of  the  electric  field  into  components  that  are  perpendicular  and  parallel  to  the 
current  flow.  This  is  to  account  for  the  surface  scattering  to  which  the  carriers  in  the  channel  are 
subjected  to.  The  resulting  Laplace's  equation,  together  with  the  "floating"  (open)  boundary 
conditions  are  then  discretised  and  solved  using  the  SLOR  method.  An  algorithm,  briefly  describing 
the  procedure,  is  shown  on  a  flow  chart  and  a  program  listing  of  the  MOSFET  split -drain  model  is 
given  in  Appendix  I.  Appendix  H  shows  the  results  of  the  two  dimensional  split -drain  simulations 
which  include  the  potential,  current  density  and  surface  charge  distributions  in  the  MOSFET 
channel  for  a  variety  of  aspect  ratios  and  drain  separations.  An  optimal  device  configuration  is 
predicted  with  its  sensitivity  studied  for  different  magnetic  field  strengths,  and  a  suitable  circuit 
that  optimises  this  sensitivity  is  shown. 

The  thesis  ends  with  Chapter  5  which  presents  a  summary  of  this  work  and  outlines  the 
possible  areas  which  require  further  investigations. 
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2.  CARRIER  TRANSPORT  IN  MAGNETIC  SEMICONDUCTORS 

When  a  semiconductor  device  is  placed  in  a  magnetic  field,  moving  charge  carriers 
(electrons  and  holes)  are  subject  to  the  Lorentz  force.  Several  phenomena  result  as  a  consequence 
of  the  Lorentz  force,  e.g.  Hall  effect,  carrier  deflection,  magnetoconcentration (MC)  and  other 
secondary  effects  not  mentioned  in  this  thesis  but  which  are  described  in  [18].  These  effects,  or  a 
combination  of  them,  may  be  exploited  in  the  design  of  transducers  for  the  detection  of  static  and 
dynamic  magnetic  fields  [2]. 

2.1  The  Current  Density  Equations 

The  presence  of  a  magnetic  field  explicitly  shows  up  in  the  transport  equations  which  now 
become,  [18], 

Jn  =  (qunnl  +  qDnVn)  -  p  n*  (?n  x  S)  ,  (2.1) 

3  =  (qa  pE  -  qDpVp)  +  u  *  (3  x  B)  .  (2.2) 

The  second  term,  the  Lorentz  force,  is  dependent  on  the  current  density,  J  and  the 

n ,  p 

->  * 
magnetic  field,  B.  The  drift  mobilities  pn  and  the  Hall  mobilities  p  n  p  are  treated  as  positive 

constants. 

Consider  just  the  electron  current  density  eqn.  (2.1).  Let  A  =  (qpn  nE  +  qDR  V  n) 
and  taking  the  vector  product  of  B  with  (2.1)  we  obtain, 

2  x  3n  =  S  x  Xn  -  pn*  (|2|)2  3n  +  un*  (S  •  3n)  S  .  (2.3) 

The  dot  product  of  B  with  (2.1)  yields 

g  .  J  =  B  •  A  .  (2.4) 

n  n 

Putting  (2.3)  and  (2.4)  into  (2.1)  and  solving  for  J  n  ,  we  arrive  at 
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f 


6 


5  = 


n 


1  +  (yn*|B|) 


'2  [K  +  M n*  (S  +  (un*)2  <S  •  Xn> 


B 


.(2.5) 


Eqn.  (2.5)  is  the  general  expression  relating  the  current  density,  Jn  and  the 
magnetic  field,  B.  We  show  two  cases  regarding  the  direction  of  the  current  density,  J  n . 
i.  When  J  n  is  perpendicular  to  B,  eqn.  (2.4)  becomes 


B 


(2.6) 


and  substituting  (2.6)  into  (2.5)  yields. 


*„x= 


1  +  (un*|B|) 


A  +  u  *  (B  x  A  ) 
n  n  n 


ii.  When  J  is  parallel  to  B,  from  eqn.  (2.1)  we  obtain. 


(2.7) 


(2.8) 


2.2  Hall  Angle 

The  Hall  angle  can  be  defined  as  the  angle  between  the  current  density  and  the  resulting 
gradient  of  the  quasi -Fermi  level. 


(a)  For  electrons. 


*  tg9Hn 


(2.9) 


where t:  is  a  unit  vector  indicating  the  direction  of  J  n  x  A  n,  0^n  is  the  Hall  angle  for  electrons, 
and  X  can  be  expressed  as  =  qpnnv(()n.  Here,  v<|>n  denotes  the  gradient  of  the 
quasi-Fermi  level  for  electrons. 
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From  (2.1)  we  get. 


3  4  = 

n  n 


n 


3  <3  •  i  )  -  S(3n  •  £  ) 

n  n  n  n 


(2.10) 


-> 


If  Jn  is  perpendicular  to  B,  using  (2.4),  B  .  =0.  Thus  from  (2.9)  and  (2.10),  we  obtain. 


*  tg0Hn  =  -“n*S  • 


(2.11) 


(b)  For  holes,  with  the  same  definition,  if  JD  is  perpendicular  to  B,  then 


*  tgeHP  ■  uP*s  • 


(2.12) 


2.3  Current  Density  Components 

•f 

If  B  is  chosen  perpendicular  to  the  device  plane  and  parallel  to  the  z-axis,  i.e. 

-y 

B  =  (0.0.B  ),  the  current  density  components  in  the  device  plane  can  be  expressed  as  the 
following  ; 


nX  1  +  <un*Bz)2 


(qy  nE  +  qD  |^) 
n  n  x  n  3x 


U  *B  (qu  nE  +  qD 
n  z  ^  n  y  ^  n  dy 


(2.13) 


^  1+^n*Bz) 


/  „  an. 

oi  nE  +  qD  — ) 
y  n  ay 


+  U  *B  (qu  nE  +  qD  — ) 
n  z  '^Mn  x  ^  n  ax 


(2.14) 


„  /  __  ,  __  an. 

=  (<*nnEz  +  qDn  3l) 


nz 


(2.15) 


PX  1  +  VV 


(qu  pE  -  qD  |£) 
^  p  x  ^  p  ax 
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J 


1 


1  +  vv2 


(2.17) 


J 


pz 


(qp  _PE  -  qD  . 

P  2  ^  p  32 


(2.18) 


Eqns.  (2.15)  and  (2.18)  are  unaffected  by  the  magnetic  field  as  expressed  by  (2.8).  A  microscopic 
treatment  of  the  transport  phenomenon  in  semiconductors  under  the  influence  of  magnetic  fields 
is  given  in  [2,19]. 

2.4  Magnetic  Field  Effects 

Consider  a  slab  of  semiconductor  that  is  uniformly  doped  n-type  silicon,  with  metallic 
contacts  at  both  ends  as  shown  in  Fig.  2.1.  Current  flowing  along  its  length  L,  in  the  y*direction, 


and  a  magnetic  field,  B  applied  in  the  z-direction  will  result  in  a  Lorentz  force  in  the  x-y  plane. 

2.4.1  Hall  effect 

If  the  geometry  of  the  slab  shown  in  Fig.  2.1  was  such  that  W/L  **  0,  i.e.L  >>  W,  a 
transverse  field  develops  across  the  slab  (x -direction)  to  compensate  for  the  Lorentz  force.  The 
transverse  field  gives  rise  to  a  Hall  voltage  [20]  set  up  across  the  slab. 

Based  on  this  effect  are  devices  realised  in  both  bipolar  and  MOS  technology.  Gallagher 
and  Corak  [21]  proposed  the  use  of  the  channel  region  of  a  MOST  as  a  Hall  plate  with  two 
electrodes  arranged  along  the  channel.  Theoretical  and  experimental  work  has  been  carried  out  on 
this  device  (see  references  contained  in  [2])  and  the  best  sensitivity  was  achieved  with  the  MOST 
operated  in  the  saturation  region. 


- 


. 
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Takamiya  and  Fujikawa  [22]  proposed  incorporating  a  differential  amplifier  and  a  Hall 
plate  in  one  device  :  the  differential  amplification  magnetic  sensor  (DAMS).  Large  sensitivities 
were  achieved  as  a  result  of  the  amplification  of  the  differential  stage;  however  the  nonlinearity 
became  large  too.  Popovic  and  Baltes  [23]  proposed  a  magnetotransistor  in  CMOS  technology,  of 
which  the  operation  is  reminiscent  of  the  SOS  magnetodiode  [24],  but  which  does  not  require  an 
AI2  O3  substrate  ;  the  high  recombining  interface  is  replaced  by  a  reverse  biased  collector 
junction  providing  a  virtually  infinite  recombination  velocity. 

2.4.2  Carrier  Deflection 

If,  on  the  other  hand,  the  geometry  of  the  slab  in  Fig.  2.1  was  such  that,  L/W  ■*  0,  i.e. 
W  >>  L,  the  transverse  field  is  absent  in  this  case,  giving  to  carrier  deflection.  Thus,  the  current 
density  vector  is  no  longer  parallel  to  the  y-axis  but  rotated  by  the  Hall  angle.  Note  that  this 
situation  is  in  fact  complementary  to  the  Hall  effect,  in  which  the  current  density  vector  is 
unaffected  and  the  electric  field  vector  is  rotated. 

Utilising  the  effect  of  carrier  deflection,  Hudson  [25]  invented  the  dual  collector 
magnetotransistor  and  used  the  difference  between  the  collector  currents  as  a  measure  of  the 
magnetic  flux  density.  The  sensitivity  figures  achieved  using  carrier  deflection  are  comparable  to  or 
as  much  as  an  order  of  magnitude  better  than  those  of  Hall  plates.  Other  advantages  of  devices 
based  on  this  principle  are  discussed  in  [25].  An  extensive  review  of  dual  and  multicollector 
magnetotransistors,  both  vertical  and  lateral  structures,  as  well  as  carrier  domain  magnetometers, 
all  based  on  the  bipolar  technology  is  given  by  [2]. 

The  principle  of  carrier  deflection  was  also  used  in  MOS  transistors  and  this  was  carried 
out  by  Fry  and  Hoey  [26].  They  used  the  Hall  electrodes  (discussed  previously  in  the  Hall  plate 
model)  as  a  pair  of  drains  additional  to  or  instead  of  the  real  drain.  Carr  and  Hong  [27]  put 
foward  a  theoretical  treatment  of  the  carrier  deflection  mechanism  in  the  dual  and  triple -drain 
MOSFET's.  Popovic  and  Baltes  presented  a  CMOS  MFS  [10]  which  consists  of  a  pair  of 
split -drain  MOS  transistors  in  a  CMOS  differential  amplifier -like  configuration.  This 
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configuration  is  reported  to  have  a  sensitivity  that  is  appreciably  higher  than  known  M05  Hall 
elements. 

2.4.3  Magnetoconcentration 

When  electrons  and  holes  are  simultaneously  injected  from  opposite  sides  of  an  intrinsic  or 
"close-to"  intrinsic  semiconductor  slab,  the  Lorentz  force  deflects  both  electrons  and  holes  to  the 
same  surface.  Since  a  Hall  voltage  is  prevented  from  building  up  (a  consequence  of  electrical 
neutrality),  the  Lorentz  force  is  compensated  by  a  carrier  concentration  gradient  perpendicular  to 
the  magnetic  and  electric  field  vectors.  The  conductivity  of  the  slab  can  be  changed  considerably  by 
allowing  recombination  of  the  carriers  at  the  surface  which  was  enhanced  (or  depleted)  of  carriers. 

Pfleiderer  [18]  presented  theoretical  and  experimental  analysis  of  the  magnetodiode  which 
is  a  device  based  on  this  effect.  Kamarinos  et  al.  [24]  proposed  a  silicon -on -sapphire  (SOS) 
magnetodiode  which  employs  a  Si-Si02  interface  on  top  for  obtaining  a  low  recombination 
velocity  (Sj^ )  and  a  Si-Al  20  3  interface  at  the  bottom  for  a  high  recombination  velocity  (s2  ). 
Crucial  requirements  to  achieve  a  high  sensitivity  are  :  a  large  surface  recombination  velocity 
difference  (s-l-s2  )  and  a  diode  thickness  that  is  much  smaller  than  the  ambipolar  diffusion 
length. 


. 


11 


Fig.  2.1  Geometrical  representation  of  a  semiconductor  slab  of 

length  L,  width  W,  and  B  is  perpendicular  to  the 


device  plane 


3.  ANALYSIS  OF  A  SEMICONDUCTOR  SLAB 


Consider  a  uniformly  doped  rectangular  silicon  slab  with  metallic  contacts  at  two  edges  as 

shown  in  Fig.  3.1.  Current  flowing  through  the  slab  in  the  y- direction  and  a  magnetic  field, 

B  =  (0,0, Bz  )  applied  perpendicular  to  the  slab  surface  will  result  in  a  Lorentz  force  in  the  x-y 

plane.  The  thickness,  h  of  the  slab  is  chosen  large  enough,  hence  any  surface  effects  at  z  =  0  and 

z  =  h  can  be  neglected.  Furthermore,  by  using  the  assumption  that  the  electrostatic  potential  ^  , 

*  * 

the  concentrations  n  and  p,  the  mobilities  pn ,  pn  ,  pp  ,  pp  are  all  independent  of  z,  we  can 
reduce  our  problem  (which  is  three  dimensional)  to  that  of  two  dimensions. 

3.1  Basic  Equations 


Similar  to  the  case  of  zero  magnetic  field,  Jnz  =  J  pz  =0.  Using  E  =  -  v  ^  and 
Einstein's  relations,  D  =  P  kT/q,  eqns.  (2.13) - (2.18)  are  now  simplified  to  two  dimensional 
current  density  equations,  viz. 


(3.1) 


(3.2) 


(3.3) 
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J  (x,  y)  = 

py 


-y 


1  +  (y  *B  ) 
P  z 


-  J  *1  (x  v)  +  kT  9 
2  t  9y  (x'  y)  +  q" 


(3.4) 


y  *B 
P  z 


1*  {Xf  y)  +  *2  i_ 

,3x  q  8x 


qp(x,  y)  | 


The  current  continuity  equations  (for  the  time  stationary  case)  and  Poisson's  equation  remain  the 
same  as  for  zero  magnetic  field  and  are  given  by  [4,5]. 


q 

•  3  =  R  , 

n 

(3.5) 

-  V 

•  3  =  R  , 

(3.6) 

q 

P 

2ip  = 

-  ^  (p-n+N) 

er 

(3.7) 

Here,  R  denotes  the  net  recombination  rate  and  N  the  net  doping  density,  -N  .  The  form  of 

eqns.  (3.1) -(3.4)  are  suitable  for  the  choice  of  the  unknown  variables  p,n  and  \p  and  will  be 
considered  as  the  fundamental  set  of  unknown  variables  throughout  the  thesis. 

In  some  cases  different  choices  of  variables  such  as  (J)^ ,  ())p  and  ip  are  used  [4,28]  in 
which  case  the  equations  (3.1) - (3.4) ,  (3.7)  must  be  suitably  modified  [4,28].  The  quasi-Fermi 
potentials  are  defined  by  [4,28], 
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exp 


kT 


U  “  <J>n) 


(3.8) 
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p  =  exp 


(3.9) 
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Fig.  3.1 


Basic  slab  geometry  :  semiconductor  slab  of  length  L  , 
width  W  ,  and  thickness,  h.  B  is  perpendicular  to 
the  drawing  plane 
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3.2  Boundary  Conditions 

In  this  section,  two  types  of  boundary  conditions  are  considered  : 


(a)  Fixed  boundary  conditions 

With  reference  to  Fig.  3.2,  the  fixed  boundaries  are  denoted  by  and  ^2^.  Along 
these  "fixed"  boundaries  the  values  p,  n  and  ^  are  given  by  infinite  contact  recombination 
velocities  and  charge  neutrality  [4,5,28].  Hence,  carriers  are  in  thermodynamic  equilibrium  (in 
terms  of  quasi -Fermi  potentials,  <l>n  =<t»p  =  Va) 


p  •  n  =  n. 


(3.10) 


p  -  n  +  N  =  0 


(3.11) 


Boundary  values  of  $  are  determined  from  the  definition  of  quasi-Fermi  potentials,  (3.8)  and 
(3.9)  and  the  space  charge  neutrality  condition  (3.11),  yielding  [4,5,28], 


•  /  N  v 

sinh  ( 2^— ) 

ni 


(3.12) 


(b)  Floating  boundary  conditions 

The  floating  boundaries  are  denoted  by  anc*  in  Fig.  3.2.  At  these  boundaries,  a 
zero  normal  component  for  the  electron  and  hole  current  densities  is  assumed  (characterising  zero 
surface  recombination)  [4,5,28]  ; 


J 

nx 
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(3.13) 
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px 


0  . 


(3.14) 
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Fig.  3.2 


Two  dimensional  slab  geometry 
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For  the  required  third  boundary  condition  with  nonzero  magnetic  field,  one  cannot  adopt  any  of 


the  conditions  commonly  used  in  the  cases  of  zero  magnetic  field.  A  zero  normal  component  of 
the  electric  field  [4,28]  would  be  inconsistent  with  the  presence  of  a  Hall  field.  If  a  zero  space 
charge  variation  normal  to  the  surface  is  asssumed  [5],  the  electron  and  hole  concentrations  at  the 
boundaries  would  be  too  tightly  coupled.  This  leads  to  the  choice  of  a  condition  that  is  somewhat 
different  than  the  conditions  commonly  used  [4,5]  along  the  floating  boundaries  ;  that  is 


3 


2 


(P 


n  +  N) 


0 


(3.15) 


However,  it  still  ensures  that  the  space  charge  is  a  smooth  function  of  the  spatial  coordinates.  A 
more  advanced  analysis  will  have  to  take  into  account  surface  states  [4,29]  and  possibly  the 
external  electric  field. 


3.3  Models  of  Physical  Parameters 

The  field  equations  contain  quantities,  e.g.,  net  recombination  rate  R,  the  mobilities  pn , 

*  * 

p  p  and  the  Hall  mobilities  p  n ,  p  ,  which  themselves  are  the  results  of  rather  complicated 
physical  mechanisms.  These  quantities  are  not  constant  but  could  depend  on  the  local  values  of 
carrier  densities,  impurity  concentrations  and  fields. 


(a)  Carrier  mobility 

In  general,  carriers  are  scattered  by  phonons  and  lattice  defects  and  in  addition  Coulomb 
scattering  at  ionized  impurity  atoms  results  in  a  further  reduction  of  the  mobility.  By  fitting 
experimental  data  [30],  an  empirical  expression  for  the  drift  mobilities  is  : 
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u  —  u 

max  mm 

1  +  <Vnref> 


(3.16) 


Here,  NT 


+  na 


is  the  total  doping  concentration,  and  the  values  of  the  constants  are 
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shown  below  : 


u  n 

UP 

NREF(cm  ) 

8.5  x  1016 

6.3  x  lo16 

a 

0.72 

0.76 

Umaxlcm2/Vs) 

1330 

495 

Wmin(crn2/Vs) 

65 

47.7 

For  the  Hall  mobility,  the  approximation 

* 

n  =  r v  (3.17) 

is  used  [19,31]  with  r  =  <x2>/<T>2,  where  <  >  denotes  the  expectation  value.  The 

parameter  x  is  the  mean  free  time  between  carrier  collisions  which  depends  on  the  carrier  energy, 
r  takes  a  value  of  1.18  for  phonon  scattering  and  1.92  in  the  case  of  ionised  impurity  scattering.  In 
all  computations  shown  in  this  thesis,  a  value  of  1.92  is  adopted  for  r. 


(b)  Recombination  and  generation 

The  bulk  recombination  rate,  R  is  assumed  to  follow  the  Shockley -Read -Hall  (SRH) 
recombination  model  [32],  viz. 


p  n  -  n. 


R  = 


n 


Ei  -  Et 

p  +  ni  exp  (— — -) 


+  T. 


E.  -  E 


(3.18) 


n  +  ni  exp ( — ^ 


Here  x  and  x  are  the  bulk  carrier  lifetimes,  E^  denotes  the  energy  level  of  trap  centers  and 
n  p  t 

E  .  the  intrinsic  Fermi  level.  Usually  E .  is  not  known  and  following  common  practice  [4,5], 
1  c 

E  =  Ei  is  assumed  in  the  computations. 
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3.4  Numerical  Method 

The  partial  differential  eqns.  (3.1) - (3.7),  (3.18),  and  the  boundary  conditions 
(3.10) -(3.15)  have  to  be  discretised  for  numerical  solutions.  In  view  of  the  rectangular  geometry 
of  the  slab,  a  finite  difference  technique  was  implemented  on  a  nonuniform,  rectangular  grid  (see 
Fig.  3.4).  This  technique  allowed  the  use  of  the  SLOR  (Succesive  Line  Over  Relaxation)  method, 
achieving  reasonably  fast  convergence  when  the  equations,  (3.1 ) - (3.7)  were  solved  simultaneously 
[5]. 

Eqns.  (3.1) -(3.7)  and  (3.18)  are  discretised  and  then  linearised  by  a  Newton  iteration 
scheme  [5].  The  linear  system  thus  obtained  are,  however,  in  many  cases  unsuitable  for  iterative 
solution  schemes  because  of  the  numerical  values  of  the  out -of  -diagonal  elements.  The  problem  of 
achieving  diagonal  dominance  can  be  removed  by  employing  the  well  known  Scharfetter-Gummel 
scheme  [33]  but  generalized  for  the  case  of  two  dimensions  [34,35]  and  nonzero  magnetic  field. 
This  generalization  is  discussed  in  the  following  section. 

3.4.1  Discretisation 

(a)  Continuity  equations 

The  analysis  of  this  section  is  carried  out  for  the  electron  continuity  equation.  The  same 
analysis  holds  for  the  continuity  equation  for  holes. 

Applying  Gauss'  theorem  to  eqn.  (3.5),  we  obtain  for  the  stationary  case, 

/  ?  •  3  dV  =  /  3  •  d3  =  /  qR  dV  .  (3.19) 

Vi  Si  Vi 

With  reference  to  Fig.  3.3,  Vj_  denotes  the  volume  element  of  the  cell  with  the  central  node  i.  The 
surface  S  ±  encloses  the  volume,  Vi  .  The  basic  assumptions  [34,35],  of  this  method  are  that  the 

-V  ->■ 

current  density  Jn  ,  electric  field,  E  and  the  diffusion  constant,  Dn  remain  constant  inside  each 
of  the  triangles  shown  in  Fig.  3.3  (e.g.  the  triangle  specified  by  the  central  node  i  and  two 
neighbouring  nodes,  j  =  1,2).  Then,  the  surface  integral  in  (3.19)  is  readily  simplified  to 
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(3.20) 


where  J  j  is  the  electron  current  density  component  flowing  from  node  i  to  a  neighbouring  node 

j’dij  is  the  length  of  the  common  boundary  between  node  i  and  node  j  (see  Fig.  3.3)  and  h 

denotes  the  thickness  of  the  slab.  To  obtain  Jn^j  eqn.  (2.7)  was  used  in  a  modified  form. 

Denoting  the  derivative  in  the  ij -direction  by  3  /  3  s,  the  derivative  normal  to  ij  (in  the 

counter-clockwise  direction)  by  3/  3a  and  the  component  of  current  normal  to  J  .  .  as  J  x.  .  . 
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(3.22) 


The  derivatives  with  respect  to  'a'  are  estimated  at  auxilliary  grid  points  A,B,C,D  (see  Fig.  3.3) 

and  these  pairs  of  points  are  denoted  by  dummy  variables  (k£)  =  (A,B),  (B,C),  (C.D),  (D,A). 

Keeping  in  mind  our  basic  assumptions  that  current  densities,  electric  field  and  mobilities  are 

constants  within  each  of  the  triangles  shown  in  Fig.  3.3,  (3.21)  and  (3.22)  are  integrated  over  s 

(between  i  and  j)  and  over  ’a’  (between  k  and  &  )  respectively.  Denoting  the  electron 

concentrations  at  i,  j,  k,Ji  as  n  ±  ,  n  ^  ,  n  k  and  n  £  and  the  electrostatic  potentials  at  i,  j,  k, l  as 

i h  ,  . ,  ij; ,  ,\b„  ,  we  obtain  from  (3.21)  and  (3.22), 

i  i  k  k 
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J  .  .  +  y  *B  J.  . 

m3  n  z  nij 
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exp 
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(3.23) 
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n 


kT  "  V 


exp[kT  <*l  '  V] 
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exp[fe  (^k  -  v]  ’  1 


(3.24) 


By  introducing  the  "Bernoulli  coefficients"  [34,35], 


B.  .  = 


13  exp 


(^i  -  q/kT 


-  ^j)  q/kT 


-1 


(3.25) 


eqns.  (3.23)  and  (3.24)  can  be  simplified  to  : 


J  .  .  +  u  *B  Jx.  .  =  j"n  .B ^ ^  -  n.B.  .1 

m3  n  z  ni]  Sij  n  [  3  J1  1  13J 


(3.26) 


r£ij  "  vn*BzJnij  =  [n*B*k  '  nkBk*] 


(3.27) 
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Eliminating  J±, .  from  (3.26)  and  (3.27),  we  obtain  for  J  . . , 

mj 


J  .  .  = 
niD 


kTy 


n 


1  +  (y  *B  ) 
n  z 


n.B.  .  -  n.B.  . 

-3  31  1  3-1  _ 

S  .  . 

ID 


y  *B 
n  z 


n£B£k  "  nkBk£ 


d.  . 

ID 


(3 


The  net  recombination  rate  (3.19)  is  assumed  to  have  an  average  value  Ri  inside  the  volume 
V  £  =  A±  x  h,  where  Ai  denotes  the  areas  of  the  square  ABCD  (see  Fig.  3.3).  Thus  eqn. 
(3.19)  becomes, 


/  J  •  dS  =  qA . hR . 

S.  n  i  1  * 

l 


(3.29) 


Substituting  (3.28)  into  (3.20)  and  equating  the  result  to  (3.29),  we  obtain. 


R.  = 
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kT 


1  +  (y  *B  ) 
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I 
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S .  .  nil  in 

D  1D  J 


(3.30) 
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Eqn.  (3.30)  shows  the  discretised  version  of  the  electron  continuity  eqn.  (3.5)  in  the  presence  of  a 
magnetic  field.  Applying  the  same  treatment  to  (3.6),  yields  the  discretised  hole  continuity 
equation  in  the  presence  of  a  magnetic  field,  viz. 
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(3.31) 
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The  advantage  of  this  discretisation  technique  is  that  it  handles  any  type  of  grid  stucture. 
For  the  rectangular  grid  structure,  shown  in  Fig.  3.3,  (3.30)  and  (3.31)  are  simplified  further  by 
linearly  interpolating  the  values  of  the  variables  p,  n  and  $  at  the  auxiliary  grid  points  A,  B,  C, 
D,  using  the  values  of  these  variables  at  the  main  grid  points.  The  interpolated  values  of  p,  n,  ^ 
and  the  simplification  of  the  Bernoulli  coefficients  are  shown  in  Appendix  A.  With  reference  to 
Fig.  3.3,  the  electron  continuity  equation  is  simplified  to  (see  Appendix  A)  ; 
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(3.32) 


The  discretised  recombination  term  ,  from  (3.18)  reads, 


p .  n . 
l 
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-  n .  . 
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p .  +  n .  . 
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(3.33) 
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In  eqns.  (3.32)  and  (3.33),  n^ 


n  3’  n4  *  ni  and  ’  ^2 


ip  ,  ip,  ,  ip.  denote  the 


electron  concentration  and  the  potential  at  the  main  points  (M,N-1),  (M  +  1,N),  (M,N  +  1), 


(M*1,N)  and  (M,N)  respectively.  Note  that  n.  denotes  the  intrinsic  carrier  concentration,  n  .  at 

ii  i 

node  (M,N).  From  Fig.  3.3,  we  note  the  following  ; 


d.,  h  (M-l)  +  h  (M) 
ii  _  x  _ x 

Sil  2hy(N-l) 

di2  hy(N-l)  +  hy(N) 
=  2hx(M) 

di3  hx(M‘1)  +  hx(M) 
Si3  "  2hy (N) 

d.  .  h  (N-l)  +  h  (N) 
i  4  y  Y 

S~ 7  =  2h  (M-l'j 

i4  x 


(3.34) 


Eqn.  (3.31)  treated  in  a  similar  fashion  yields  the  hole  continuity  equation  as  ; 
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(3.35) 


where  p  ,  p2  ,  P3  ,  P4  and  p  .  denote  the  hole  concentrations  at  the  main  points  (M,N-1), 
(M  +  l.N),  (M.N  +  l),  (M-l.N)  and  (M.N)  respectively.  The  term,  R.  is  given  by  eqn.  (3.33). 
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Fig.  3.3  A  non-uniform  rectangular  grid 
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(b)  Poisson's  equation 

Using  a  finite  difference  technique,  eqn.  (3.7)  is  discretised  by  estimating  the  variables  p.n 
and  ip  at  the  main  points  (see  Fig.  3.3),  and  we  obtain  as  a  direct  consequence  of  eqn.  (3.7),  the 
discretised  version  of  Poisson's  equation  : 


tiy(N-l)  |hy(N)  +  h  (N-l)  j  *1  +  n  (M)  [h  (M)  +  h  (M-l)]  *2 

••  A  *  — 


hy(N)  Lhy  (N)  +  hy  (N-l)  J  *3  hx(M)  [hx(M)  +  hx(M-l)J  *4 


h  (M)  h  (M-l) 

A  A 


h  (N)  h  (N-l) 


(3.36) 


+  -3—  p . 

cSi  1 


-2-  n. 
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=  _  _SL 


N 


£Si  1 


where  N  .  denotes  the  net  doping  at  node  (M,N) , 

1 


(c)  Floating  boundary  conditions 

Eqns.  (3.13) - (3.15)  are  discretised  on  the  nonuniform  grid  shown  in  Fig.  3.4,  for  the 
boundaries  ©  and  ©  (see  Fig.  3.2).  At  these  floating  boundaries,  the  variables  p,  n  and 
at  the  meshpoint  on  the  floating  line  are  expressed  in  terms  of  the  variables  at  the  neighbouring 
inner  points  [5].  This  is  described  in  the  following  section  : 

From  eqn.  (2.13),  the  discretised  version  of  (3.13)  for  boundary  (see  Fig.  3.2)  reads  ; 


n(l,  N)  =  n (2,  N)  exp  (  3_  [*(1,  N)  -  *(2,  N)  ]  +  kR  }  , 


(3.37) 


■ 
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where. 


N-l)  -  ip(2,  N+l)  ] 


In 


n  ( 2 ,  N-l)  ) 
n  (  2 ,  N+l)  f 


(3.38) 


and  b 

n 


un*Bzhx(1)/thy(N)  +  hy<N-l)] 


(3.39) 


From  eqn.  (2.16),  the  discretised  version  of  (3.14)  reads  ; 

p(l,  N)  =  p  ( 2 ,  N)  exp|^|  [i)i  ( !,  N)  -  ip  (2,  N)  ]  +  kp  1  ,  (3.40) 


where, 

kp  =  bp{fe  (3.41) 


and 

bp  =  up*Bzhx(l)/[hy(N)  +  hy(N-l)]  .  (3.42) 


The  discretised  version  of  (3.15)  reads  ; 
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As  eqns.  (3.37)  and  (3.40)  are  nonlinear  in  p,  n  and  \p  ,  we  express 

p(M,N)  =  p^M.N)  +  5p(M,N),  n(M,N)  =  iftM.N)  +  6n(M,N), 

(M.N)  =  ^M,N)  +  5ip(M,N)  and  apply  Taylor's  expansion  theorem  [5]  to  the  nonlinear 


m 
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functions,  neglecting  higher  order  terms.  The  linearised  versions  of  (3.37)  and  (3.40)  are  put  into 
(3.43),  hence  obtaining  the  linearised  'P  (l.N).  We  realise  that  following  this  procedure,  the 
linearisation  is  introduced  prior  to  applying  the  boundary  conditions.  In  general,  this  may  lead  to 
errors  due  to  subtractive  cancellation;  however,  we  have  verified  that  the  results  obtained  by 
Kurata's  procedure  [5]  agree  with  those  obtained  by  applying  the  boundary  conditions  prior  to 
linearisation. 

The  linearised  versions  of  eqns.  (3.37)  and  (3.40),  together  with  eqn.  (3.43)  can  be 
summarised  in  matrix  vector  notation  as 


9(1,  N)  =  Gt(2,  N)  6(2,  N+l)  +  HT ( 2 ,  N)  9(2,  N-l) 

+  It(2,  N)  6(2,  N)  +  JT(2,  N)  6(3,  N)  +  K(N) 


(3.44) 


where  in  general, 


0(X,  N)  =  [p  (X,  N)  n(X,  N)  ip  (X,  N)  ] 


(3.45) 


and  GT  ,  HT  ,  IT  ,  JT  are  3  x  3  matrices  whose  elements  are  given  in  Appendix  B.  The  constant 
vector  K(N)  is  a  1  x  3  matrix  whose  elements  are  given  in  Appendix  B. 

In  the  same  way,  we  obtain  a  matrix  vector  equation  for  boundary  (T ^  as  (see  Fig.  3.2) 


e(XMAX,  N)  =  Gt(XMAX-1,  N)  6 (XMAX-1,  N+l) 


+  HT (XMAX-1,  N)  6 (XMAX-1,  N-l) 

+  It(XMAX-1,  N)  6 (XMAX-1,  N) 

+  JT (XMAX-1,  N)  6(XMAX-2,  N)  +  L (N)  . 

The  elements  of  the  matrices,  GT  ,  HT  ,  I  T  ,  JT  and  L  for  this  floating  boundary 


(3.46) 


are  in 


Appendix  B. 
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Fig.  3.4  Discretisation  of  equations  (3. 13) -(3. 15)  at  the  floating 

boundaries  (7)  and  (7) 
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3.4.2  Linearisation  of  the  Field  Equations 

As  eqns.  (3.5) -(3.7)  are  nonlinear  with  respect  to  the  fundamental  variables  p,  n  and  ip  , 
Newton  s  method  is  used  to  linearise  these  equations,  leading  to  the  simultaneous  solution  of  p,  n 
and  \p  [5].  As  in  [4,35],  the  eqns.  (3.5)-(3.7)  can  be  written  in  the  form: 

Fp(Pr  n,  1)1 )  =  i  V  •  5  +  R  =  0  f  (3.47) 

Fn(p,  n,  +)  =  -  |  V  •  3n  +  R  =  0  ,  (3.48) 

9  rr 

F ,  (P/  n,  ip)  =  v  ip  +  -2—  (p-n+N)  =  0  ,  (3.49) 

^  Si 


where  the  RHS  of  eqns.  (3. 47) -(3. 49)  are  given  by  (3. 32) -(3. 36).  We  can  express  p,  n  and  ^  as 
p  =  +  6p  ,  n  =  +  6n  ,  ip  =  ip  +  S\p  , 

where,  p °,  n°  and  \p°  are  trial  values  (or  values  from  a  previous  iteration)  and  6p,  6 n, 
represent  the  'correction'  or  'improvement'  on  p°  ,  n°  ,  \p  °  respectively,  through  a  Newton 
iteration  step.  Applying  Taylor's  expansion  theorem  to  the  nonlinear  functions,  neglecting  the 
higher  order  terms,  we  get 
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Here, 


3  / 

(Pr  n 


n,  if;) 

¥) 


is  the  Jacobian  matrix  evaluated  at  (p,  n,^)k  and  (  6p,  6n,  6i|;)k+1  is  the  correction  at 
iteration  step  (k  +  1),  given  by 


Pk+1 


k  k+1 

p  +  6p 


/ 


nk+1 


6  n 


k+1 


(3.51) 


k+1 


+  6  ip 


k+1 


/ 


1c  Ic  1c 

where  k  =  0,  1,  2...  and  p  ,  n  ,  ip  are  the  initial  guess  values.  Eqn.  (3.50)  represents  the 
basic  equation  for  a  single  node  (meshpoint).  Taking  into  account  all  the  five  nodes  (see  Fig.  3.3), 
where  the  coefficients  and  constants  carry  different  values  at  each  of  the  nodes,  we  can  summarise 
the  two  dimensional  problem  to  yield  a  matrix  vector  equation  that  treats  p,  n  and  i p  as  our 
unknowns  [5]  ; 


AT(M,  n)  e(M,  N-1)  +  BT(M,  n)  e(M,  n)  +  ct(m,  n)  e (m,  n+1) 

(3.52) 

+  D  (M,  N)  0 (M-l,  N)  +  ET(M,  N)  0(M+1,  N)  =  FT (M,  N)  , 


where 

6  (M,  N)  =  0°(M,  N)  +  60  (M,  N) 


P(M,  N) 
n(M,  N) 

<l>  (M,  N) 


✓ 
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The  definitions  of  the  matrices  and  vectors ; 


A^M,  N)  =  a  , 

bt(m, 

N)  =  b..  , 

cT(M, 

N)  =  c. .  , 

ID 

dt(m,  N)  =  ai;j  , 

et(m. 

N)  =  e.  , 

Ft(M, 

N)  =  fi  (3.53) 

i  “  1  r  3  r  j—  1/3  / 

whose  elements  are  related  to  the  coefficients  of  the  continuity  and  Poisson's  equations,  are  given 
in  Appendix  C. 

3.4.3  Solution  Procedure 

Eqn.  (3.52)  is  solved  for  every  meshpoint  in  our  nonuniform  grid  with  the  necessary 
modifications  by  (3.45)  and  (3.46)  at  the  floating  boundaries.  For  example,  with  a  grid  size  of 
1  £  M  1  50,  1  g  N  g  50,  we  have  2500  equations  that  must  be  solved  simultaneously.  Since 
there  are  three  unknown  variables  at  each  meshpoint,  the  total  number  of  unknown  variables  for 
our  two  dimensional  system  is  7500.  Direct  methods  such  as  Gaussian  elimination  to  solve  such  a 
system  will  not  be  appropriate.  This  led  to  the  choice  of  the  SLOR  (Succesive  Line 
Over  Relaxation)  method  [5,37]  to  solve  the  linear  system  of  equations  (3.52).  The  SLOR  is  an 
iterative  method,  one  in  which  the  computation  must  be  iterated  until  an  accurate  solution  for  the 
system  of  linear  equations  (3.52)  is  obtained. 

Once  the  linear  system  of  equations  is  solved,  a  single  Newton  iteration  is  over.  Another 
Newton  step  is  earned  out  where  the  matrix  vector  coefficients  of  eqns.  (3.52)  is  updated  with  the 
most  recent  values  of  the  variables.  The  Newton  iterations  are  continued  until  the  original  system 
of  nonlinear  equations  is  solved  sufficiently  well.  The  Newton -SLOR  method  described  consists  of 
outer  iterations  (Newton)  and  inner  iterations  (SLOR).  Fig.  3.5  illustrates  a  flow  chart  description 


of  the  method. 
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A  specific  feature  of  the  SLOR  method  is  to  regard  the  unknown  quantities  on  two 
neighbouring  points  0(M,N-1)  and  0(M,N  +  1)  as  already  known  quantities  [5].  This  simplifies 
eqn.  (3.52)  to 

Bt(M,  N)  e(M,  N)  +  Dt(M,  N)  0 (M-l ,  N)  +  E  (M,  N)  0(M+1,  N) 

=  F(M)  for  1  ^  M  ^  XMAX-1  ,  (3.54) 

where, 

F(M)  =  Ft(M,  N)  -  At(M,  N)  0 (M,  N-l)  -  C  (M,  N)  0  (M,  N+l)  . 

Eqn.  (3.54)  is  solved  with  the  necessary  modification  in  coefficients  [5]  according  to  the  boundary 

conditions  (3. 44) *(3.46),  to  fit  each  individual  case. 

Using  a  nonuniform  grid  of  50  x  50  nodes  and  working  with  double  precision  arithmetic, 

the  computation  of  one  solution  required  approximately  2Mbytes  of  memory  and  50-200  secs,  of 

CPU  time  on  the  Amdahl  580  computer.  For  every  Newton  cycle,  15-30  SLOR  iterations  were 

-4 

implemented.  For  a  relative  precision  of  10  typically  8-20  Newton  cycles  were  required.  In  the 
SLOR  iterations,  the  relaxation  parameter,  w  (see  Fig.  3.5)  varied  between  0.1  and  1.9  depending 
on  the  device  parameters  and  the  initial  guess  values  of  0.  Depending  on  the  initial  guess,  the 
computations  may  require  up  to  one  third  more  CPU  time  than  the  case  of  zero  magnetic  field 
[38]. 

To  check  if  the  true  solution  was  actually  obtained  for  the  given  equations,  the  following 
tests  were  made  [5] ; 

i.  both  sides  of  eqn.  (3.52)  were  checked  for  agreement  at  desired  meshpoints  to  see  if  the  line 
equations  are  really  solved  through  the  SLOR  method. 

ii.  the  discretised  nonlinear  eqns.  (3.32),  (3.35)  and  (3.36)  were  checked  at  desired  meshpoints 
to  see  if  there  was  agreement  between  the  left  and  right  hand  sides. 

iii.  check  for  conservation  of  the  terminal  currents. 
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Fig.  3.5 


Flow-chart  of  the  Newton -SLOR  method 
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3.5  Results  and  Discussion 

In  this  section,  distributions  of  interest,  such  as  equipotential  lines,  current  lines  and 
carrier  concentrations  are  shown  for  uniformly  doped  rectangular  slabs  of  different  slab 
geometries  and  for  different  magnetic  field  strengths,  the  magnetic  field  being  perpendicular  to  the 
device  plane. 

The  electric  potential  and  current  densities  were  computed  for  a  square  silicon  slab  with 
metal  contacts  and  with  doping  levels  of  10  16  cm  for  both  p  and  n*type  material  [16].  Plots 
of  the  equipotential  curves  and  current  lines  in  the  slab  are  shown  for  magnetic  field  strengths  of 
8  =  1  t  and  10  t  (see  Figs.  3.6  *  3.9).  Equipotential  lines  are  parallel  close  to  the  metal  boundary 
and  the  current  lines  are  parallel  close  to  the  'floating'  (open)  boundary.  But  far  from  the 
boundaries,  however,  the  physical  mechanisms  in  the  slab  cannot  be  visualised  in  terms  of  Lorentz 
deflection  nor  Hall  voltage  ;  both  these  mechanisms  are  involved  in  a  complex  way  [16]. 

The  dependence  of  the  electric  potential  and  current  density  on  the  device  geometry  is 

shown  in  [17].  Computations  were  made  for  various  geometries  and  plots  of  the  equipotential 

curves  and  current  lines  are  shown  for  two  cases,  W/L  =  1/4  and  W/L  =  4  (see  Figs.  3.10  and 

3.11).  In  the  case  of  the  long  slab  (W/L  =  1/4),  current  lines  are  parallel  to  the  open  slab 

boundaries  except  for  regions  close  to  the  metal  contacts.  There  is  a  constant  Hall  field,  E  x  across 

the  width  of  the  slab.  In  the  case  of  the  short  slab  (with  wide  contacts,  W/L  =  4),  carrier 

deflection  prevails  in  most  of  the  device  area  [17].  There  is  no  presence  of  the  Hall  field  (E  x  =  0) 

and  the  current  lines  are  skewed  by  the  Hall  angle.  The  variations  of  the  electric  field  component 

ratio  Ex  /E  and  the  current  density  components  Jnx  /Jny  across  the  width  of  the  slab  for  the 

geometries  discussed  is  shown  in  Figs.  3.12  and  3.13.  Both  these  ratios  have  a  maximum  value  of 
£ 

tan  ©  Hn  =  p  n  Bz  *  where  tan  0Hn  is  the  Hal1  angle  for  electrons* 

In  our  discussion  so  far,  we  have  only  considered  doped  material,  i.e.  when  the  doping 

level  is  much  greater  than  the  intrinsic  concentration.  Under  such  conditions  [16,17],  the 
magnetoconcentration  effects  were  insignificant.  However,  if  the  doping  level  is  much  smaller  than 
the  intrinsic  concentration  [38],  the  magnetoconcentration  dominates  the  hole  and  electron 
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densities  as  shown  in  Figs.  3.14  (a)  and  3.14  (b).  The  intrinsic  condition  is  modeled  by  an 

o 

assumed  elevated  temperature  of  500  K  using  the  expression  [5] 
ni(T)  =  1.22  x  10  16  T  3//2exp(-1.15  /  2kT).  At  500  K,  the  intrinsic  concentration  is 
approximately  2.2  x  10 /  cm J .  Fig.  3.14  (c)  shows  the  combined  effects  of 
magnetoconcentration  and  metal  boundaries  resulting  in  concentration  gradients  in  both  x  and 
y- directions,  yielding  an  unusual  space  charge  distribution.  Under  such  conditions  [38],  the 
equipotential  lines  show  only  small  gradients  (i.e.  small  Hall  field  components),  while  the  current 
lines  crowd  on  the  high  concentration  side  of  the  device  indicating  a  local  conductivity  change  (see 
Fig.  3.14  (d». 
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Fig.  3.6  Equipotential  curves  for  n-type  Si  device, 

V  =  0.1  V,  L  =  W  =  100  pm  and  B  =  1  t 
a 


Fig.  3.7 


Current  lines  corresponding  to  Fig.  3.6 
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Fig.  3.8 


Equipotential 

V  =  0.1  V,  L 
a 


curves  for  p-type  Si  device. 


=  W  =  100  pm  and  B  =  10  t 


Fig.  3.9 


Current 


lines  corresponding  to  Fig.  3.8 
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Fig.  3.10  (a)  Equipotential  and  current  lines  for  an 

n-type  Si  device,  L  =  4W,  B  =  1  t,  =  1  V 


(b)  Same  as  (a)  but  -  2  t 


0 


Fig. 


3.11  (a)  Equipotential  and  current  lines  for 

n— type  Si  device,  L  =  W/4,  B  =  1  t 

z 

(b)  Equipotential  and  current  lines  for 

Si  device,  L  =  W/4,  B  =  2  t 

z 


an 


a 


p-type 


'  ■ 
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Fig.  3.12  Electric  field  component  ratio  E^/E^  across 

the  middle  (y  =  L/2)  of  the  devices  corresponding 
to  Figs.  3.6,  3.10  (a)  and  3.11  (a) 


Fig.  3.13  Current  deflection  Jx/J  across  the  middle 

(y  =  L/2)  of  the  devices  corresponding  to  Fig. 


3.12 
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Fig.  3.14  Intrinsic  Si  square  device  ; 

L  =  W  =  4  ym,  V  =  1  V,  B  =  1  t,  T  =  500  K 

cl  Z 

(a)  Hole  concentration 

(b)  Electron  concentration 

(c)  Space  charge  distribution 

(d)  Equipotential  and  current  lines 
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4.  ANALYSIS  OF  THE  MOS  MAGNETIC  FIELD  SENSOR 
Consider  an  N- channel  MOSFET  operated  in  the  linear  region  as  shown  in  Fig.  4.1.  The 
magnetic  field,  B  is  assumed  to  be  perpendicular  to  the  channel,  i.e.  B  =  (0,0, B7  )  ;  thus  the 
Lorentz  force  acts  in  the  plane  of  the  channel  (x-y  plane).  As  in  the  case  of  zero  magnetic  field, 
we  assume  no  current  flow  in  the  z- direction. 

4.1  Assumptions 

In  the  following,  the  assumptions  made  in  the  analysis  of  the  MOSFET  are  discussed 

[4.31]  : 

(a)  the  gate  structure  corresponds  to  an  ideal  MOS  diode  so  that  any  interface  traps,  fixed  oxide 
charges  or  work  function  differences  can  be  neglected.  In  a  nonideal  gate,  the  effect  of  fixed 
charges  and  difference  in  work  functions  cause  a  voltage  shift  corresponding  to  the  flatband 
voltage,  V  FB  .  This  in  turn  causes  a  change  in  the  threshold  voltage,  VT  . 

(b)  the  doping  in  the  channel  is  uniform. 

(c)  the  drift  and  Hall  mobilities  in  the  inversion  layer  are  constants  for  given  gate  and  drain  bias 

9y  9y* 

voltages.  This  is  done  so  by  neglecting  — n ,  — n  across  the  channel  and  assuming  no  lateral 

3z  3z 

dependence  of  y  ,  y  n  (for  MOSFET's  operating  in  the  linear  region  under  low  gate  bias 

voltages,  E  (x,y)  can  be  assumed  to  be  constant), 
z 

(d)  the  reverse  leakage  current  is  negligibly  small. 

(e)  the  transverse  field,  E  »  E  ,  the  longitudinal  field  in  the  channel,  corresponding  to  the 

z  y 

so-called  Gradual  Channel  Approximation.  This  means  that  the  charges  in  the  system  are  governed 
by  the  vertical  field  only. 

(f )  within  the  inversion  layer,  the  hole  current  density  can  be  neglected  since  f  «  t  n  for 

V  >0  and  V  <0.  The  net  recombination  rate  can  be  assumed  zero  in  the  channel.  This 
D  —  B 

yields  the  drain  current,  I  to  be  an  electron  current  of  vanishing  divergence. 
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OXIDE 


METAL 

ELECTRODE 


Fig.  4.1 


Schematic  of  an  N-channel  MOSFET 
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4.2  Basic  Equations 

Following  these  assumptions  and  considering  only  drift  currents,  the  charge  in  the 
inversion  layer  is  given  by  [31] 


Qn (x'  y)  =  Cox  1V(X'  +  2*b  '  VG^ 

+  ^si^A  Lv(Xf  y)  +  2*bJ 


(4.1) 


Here,  V  (x,v)  denotes  the  reverse  bias  voltage  between  a  point  (x,y)  in  the  channel  and  the  source 
electrode  (which  is  assumed  grounded).  C  is  the  gate  oxide  capacitance  per  unit  area  ; 


C 


ox 


£Si02/t 


ox 


(4.2) 


where  t  is  the  oxide  thickness.  $  is  the  potential  difference  between  the  bulk  Fermi  level 
ox  B 

(Epp  )  and  the  intrinsic  Fermi  level  (E  i ).  \p  is  given  by 


log  e 


(4.3) 


where  kT/q  =  25.7  mV  at  room  temperature,  NA  the  bulk  doping  density  and  np  the  intrinsic 
concentration.  As  an  approximation,  the  surface  potential  \p  q  at  inversion  can  be  expressed  as 


[31] 


4>s  (x,  y)  -  2\Pb  +  V  (x,  y) 


(4.4) 


In  the  following,  we  discuss  a  two  dimensional  approximation  to  the  exact  solution  which 

is  obtained  by  integrating  the  electron  current  density  (3.1)  and  (3.2)  over  the  channel  depth,  t. 

This  was  carried  out  using  the  assumptions  discussed  previously  (in  particular  assumption  (c))  and 

3  V  3  V 

further  neglecting  the  dependence  of  Ty  on  z-  The  integration  yields  two  dimensional  line 
current  densities  (current  per  unit  length)  and  they  read  ; 


' 
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j  (x,  y)  = 


-l 


nx 


1  +  K /BJ 

n  z 


3V  ,  ,  3  V, 

y„  t—  +  u  *B  (y  t— ) 
n  3x  n  z  n  3y 


Qn(x,  y)f(4.5) 


Jny <x'  y)  = 


-1 


1  +  (y  *B  ) 
n  z 


3V  *n  i  9Vi 

y  —  —  y  *B  (y  — — ) 
n  8y  Mn  z  n  2x 


Qn  (x,  y).  (4.6) 


Here, 

~  t 

Jnx(x'  y)  “  /  Jnx(x'  y'  z)  dz  ' 

t 

Jny(x'  y)  =  J  Jny(x'  y'  z)  dz  ' 

t 

Qn(x'  y)  =  q  /  n  (x,  y,  z)  dz  . 


(4.7) 


(4.8) 


(4.9) 


The  integration  is  carried  out  from  the  semiconductor -oxide  interface  (z  =  0)  to  the  point  (z  =  t) 
at  which  E  i  intersects  E Fn  (the  quasi -Fermi  level  for  electrons). 

Integrating  (4.1)  from  the  source  (0  V)  to  a  potential  V,  at  any  arbitrary  point  (x,y)  in 
the  channel  [39],  viz. 

Y  i 

J  Qn(V' )  dV’  =  Qn(V)  ,  (4.10) 


where  V'  is  a  dummy  variable.  From  (4.10)  we  note  the  following  ; 


Q  (V)  = 
3x  n 


3$  (V) 
n 

3x 


3V 

3y 


Qn(v) 


3&n(V) 

”9y 


(4.11) 


(4.12) 


- 


■  ,  '  ,T 
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Substituting  (4.11)  and  (4.12)  into  (4.5)  and  (4.6)  and  introducing  a  new  function 


u  =  f  Qn<V)  -  an(V) 


(4.13) 


we  can  simplify  (4.5),  (4.6)  as  follows  : 


J  (x,  y) 
nx  '  x ' 


n 


1  +  (u  *B  ) 
n  z 


3ll 

- —  —  u  *B  - 

3x  n  z  By 


(4.14) 


J  (x,  y)  = 
ny  J 


u 


n 


1  +  (u  *B  ) 
n  z 


■|^  +  y  *B  — 
By  n  z  3x 


(4.15) 


The  electron  continuity  eqn.  (3.5),  using  assumption  (f)  becomes 


J  (x,  y)  +  J  (x,  y)  =  0 

3  x  nx  J  By  ny  J 


(4.16) 


Hence  from  (4.14),  (4.15)  and  (4.16)  we  arrive  at 


2  2 
8 u 


=  0 


3x‘ 


By' 


(4.17) 


Note  that  to  arrive  at  (4.17),  we  assumed  that  diffusion  current  is  negligible  along  the  channel,  u 
can  be  determined  from  (4.1),  (4.10)  and  (4.13)  to  read 


u  =  f  <2*B  -  V  Cox  +  Cox  <1?  -  2*B  +  V  V 


¥  V2  +  IF  (2cSiNAqV  +  4eSiNAq*B)1/2 


<2esiNAqV  +  4^*^) 

b  1  A 


3/2 


(4.18) 


•• 
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4.3  Boundary  Conditions 

Eqn.  (4.17),  A  u  =  0  requires  boundary  conditions  in  order  to  completely  specify  the 
boundary  value  problem  to  be  solved.  With  reference  to  Fig.  4.2  the  boundary  conditions  are  the 
following  : 


(a)  Ohmic  contacts 


At  these  contacts  (drain  and  source  electrodes)  infinite  contact  recombination  velocity  (3.10)  and 
space  charge  neutrality  (3.11)  are  assumed,  hence  carriers  are  in  thermodynamic  equilibrium.  At 
the  drains,  eqn.  (4.18)  yields  : 


u  ,  =  u 
Dl 


V  =  V 


(4.19) 


Dl 


At  the  source  (which  is  grounded),  using  (4.18)  we  obtain  ; 


0  V 


(4.20) 


(4.21) 


(b)  Floating  boundaries 

At  these  boundaries  (denoted  by  Q,  Q  and  0 
electron  current  density  is  assumed  to  be  zero,  yielding 


in  Fig.  4.2)  the  normal  component  of 


(4.22) 


For  the  region  between  the  drains,  i.e.  the  floating  boundary  denoted  by 
we  apply  condition  (4.22),  Jny.  =0,  which  from  (4.15)  yields 


in  Fig.  4.2, 
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9u 

9y 


u 


9u 

9x 


(4.23) 


Along  the  remaining  floating  boundaries 


condition  (4.22),  J  ^  =0,  which  from  (4.14)  yields 


shown  in  Fig.  4.2,  we  maintain  the 


9u  _  9u 

—  u  *B  - 

9x  n  z  9y 


(4.24) 


Fig.  4.2 


MOSFET  channel  geometry  for  a  split -drain 
configuration 
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4.4  Carrier  Mobility 

Apart  from  the  usual  scattering  mechanisms,  mentioned  in  the  previous  chapter,  carriers 
in  the  channel  of  MOSFET  devices  are  further  subject  to  surface  scattering  or  surface  roughness 
which  strongly  degrades  the  mobility  even  further.  Following  [40]  this  effect  is  taken  into  account 
by  assuming  the  mobility  to  be  a  product  of  two  functions,  each  of  which  includes  a  field 
component  perpendicular  or  parallel  to  the  current  flow  [40],  viz. 

Un(NA'  Ey'  V  =  v  0  f(NA'  V  g(Ez)  *  (4.25) 


Here,  Ez  is  the  field  component  perpendicular  to  the  current  flow  and  Ey  is  the  component  of 

field  parallel  to  the  current  flow'.  pQ  denotes  the  low  field  mobility  of  electrons  in  intrinsic  Si  and 

2 

a  value  of  1500  cm  /  Vs  is  assumed  [31]. 

The  function  f(NA,Ey)  represents  the  property  of  bulk  mobility  in  Si  and  is 
approximately  given  by  [33], 


f  (N 


A' 


E  )  = 

y 


N. 


l  + 


Vs 


+  N 


(E  /A)  ' 

yA  + 


E 


2  -1-1/2 


(4.26) 


Here,  S,  N,  A,  F  and  B  are  constants  whose  values  are  given  below. 


s 

N 

A 

F 

B 

350 

_  , n  16  /  3 

3  x  10  /cm 

3.5  x  10  3  V  /  cm 

8.8 

7.4  x  10  3  V  /  cm 

Na  is  the  bulk  doping  density  and  a  constant  Ey  field  is  assumed  in  the  channel,  viz. 


E  = 

y 


<vD  -  Vs)/L 


(4.27) 


. 
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where  L  is  the  channel  length. 

The  function  g(Ez)  represents  the  property  of  surface  mobility  [40].  Surface  mobility 
decreases  with  increasing  gate  field  and  this  property  can  be  approximately  expressed  as  [40] 


g  =  (1  +  aE  )~  ly/2 
z 


(4.28) 


For  an  N -channel  M05FET,  a  =  1.539  x  10  ^cm  /  V  and  E  is  averaeed  over  the  channel 

z 

length  to  obtain 


E 

z 


0  V 

2*b  "  2"> 


(4.29) 


For  the  Hall  mobility. 


the  approximation  (3.17)  is  used. 


4.5  Numerical  Method 

The  partial  differential  eqn.  (4.17),  A  u  =  0  and  the  boundary  conditions  (4.23)  and 
(4.24)  are  descretised  and  numerically  solved  using  the  SLOR  technique  (discussed  in  the  previous 
chapter).  As  the  geometry  of  the  MCJSFET  channel  is  rectangular,  the  discretisation  was  based  on 
a  finite  difference  scheme  implemented  on  a  nonuniform  grid  (as  shown  in  Fig.  3.3). 


4.5.1  Discretisation 

(a)  Laplace's  Equation 

Following  the  same  technique  as  done  in  obtaining  the  discretised  Poisson's  eqn.  (3.36), 


the  discretised  A  u  =  0  reads  ; 


' 


* 
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2 

hy(N-l)  thy(N)  +  h  (N-l)]  u(M'  N_1) 

_  2  2  1 

h  (M) h  (M-l)  +  h  (N)  +  h  (N-l)  u(M'  N) 
x  y  y 

,  _  2 

hy (N) [hy (N)  +  h  (N-l)]  N+l)  (4.30) 

2 

h  (M-l) [h  (M)  +  h'TM-1) I  u(M"1'  N) 

x  *■  X  X  J 

_ 2 

+  E  (M)  [h  (Ml  +  h  (M-l)  ]  u(^4+lf  N)  -  0 

XX  X 


Eqn.  (4.30)  yields  a  vector  equation  for  the  two  dimensional  problem  : 

At(M,  N)  u(M,  N-l)  +  Bt(M,  N)  u(M,  N)  +  Ct(M,  N)  u (M,  N+l) 

(4.31) 

+  Dt(M,  N)  u (M-l ,  N)  +  ET(Mr  N)  u(M+l,  N)  =  FT (M,  N) 

with  definitions  of  vectors  given  in  Appendix  G. 


(b)  Boundary  Conditions 

At  the  fixed  boundaries,  the  value  of  u  is  definite  so  that  the  relevant  terms  in  eqn.  (4.31) 
are  transferred  to  the  right-hand  side  to  be  added  to  the  constant  term  FT  (M.N). 

At  the  floating  boundaries,  the  variable  u  at  a  meshpoint  on  the  floating  line  is  expressed 
by  those  at  the  neighbouring  inner  points  [5].  This  is  described  in  the  following  section: 

(i)  At  boundaries  ©  and  (^2^.  we  discretise  eqn.  (4.24),  (see  Figs.  4.3  and  4.4)  and 
summarising  the  results  in  vector  form,  as  : 


/ 


■ 


■ 


' 


u(l,  N)  Gt(2,  N)  u(2,  N+l)  +  Ht(2,  N)  u(2,  N-l) 
+  It(2,  N)  u  ( 2 ,  N)  +  Jt(2,  N)  u ( 3,  N)  , 


54 


(4.32) 


drain 

LJ 

till  1 

Y 

(1,NV 

(2 ,  N ) 

(3,  N) 

float  Ing 
boundary 

© 

(2,  N-1 ) 

LJ - 

1-J-l  J... 

J  l  l 

source 


Fig.  4.3 


Discretisation  of  (4.24)  at  (T) 
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u (XMAX,  N)  = 

Gt(XMAX-1, 

N) 

u (XMAX-1 , 

N+l) 

+ 

Ht(XMAX-1, 

N) 

u (XMAX-1, 

N-l) 

+ 

It(XMAX-1, 

N) 

u (XMAX-1, 

N) 

+ 

Jt(XMAX-1, 

N) 

u (XMAX-2 , 

N) 

(4.33) 


The  coefficients  of  the  variables  at  the  inner  mesh  points  are  given  in  Appendix  G. 


Y 
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/  /  /  /  /  /  /  / 
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© 


source 


Fig.  4.4 


Discretisation  of  (4.24)  at  (T) 
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(ii)  At  the  region  between  the  drains  (denoted  by  in  Fig.  4.2)  we  discretise  eqn.  (4.23),  (see 
Fig.  4.5)  and  summarising  the  results  in  vector  form  : 


u(M,  YMAX)  =  P  (M,  YMAX-1)  u(M+l,  YMAX-1) 
+  QT(M,  YMAX-1)  u (M-l,  YMAX-1) 
+  R  (M,  YMAX-1)  u (M,  YMAX-1) 

+  S  (M,  YMAX-1)  u (M,  YMAX- 2 ) 


(4.34) 


The  coefficients  of  the  variables  at  the  inner  mesh  points,  i.e.  P  ,  Q_  ,  R  and  S  are  given  in 

XT  T  T 

Appendix  G.  The  variables  at  the  mesh  points  on  the  floating  boundary  lines  (specified  by  regions 
©•©  and  are  eliminated  from  eqn.  (4.31)  by  (4.32),  (4.33)  and  (4.34). 
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Discretisation  of  (4.23)  at  (T) 
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4.5.2  Solution  Procedure 

The  solution  of  eqn.(4.31)  is  a  very  large  problem  because  of  the  number  of  mesh  points 
in  both  the  x  and  y-directions.  e.g.  1  =  M  5  50  and  1  =  N  =50  would  result  in  2500  equations 
being  solved  simultaneously.  Hence,  the  choice  of  an  iterative  method  (SLOR)  [5,37],  as  opposed 
to  a  direct  solution  method  such  as  Gaussian  elimination. 

A  specific  feature  of  the  SLOR  method  is  to  regard  the  unknown  quantities  on  two 
neighbouring  points  u(M,N-l)  and  u(M,N  +  l)  as  already  known.  Hence,  eqn.  (4.31)  becomes, 

Bt(M,  N)  u(M,  N)  +  D  (M,  N)  u(M-l,  N) 

(4.35) 

+  Et(M,  N)  u(M+l,  N)  =  F(M)  , 

where 

F (M)  =  Ft(M,  N)  -  At(M,  N)  u (M,  N-l)  -  CT(M,  N)  u(M,  N+l) 
for  1  ^  M  <  XMAX-1  . 

Eqn.  (4.35)  is  solved  with  the  necessary  modification  in  coefficients  [5]  according  to  the  boundary 
conditions  to  fit  each  individual  case.  The  execution  of  the  SLOR  method  is  described  in  the  flow 
chart  shown  in  Fig.  4.6. 

The  actual  computation  time  and  the  number  of  SLOR  cycles  depend  on  the  initial  guess 

values  of  u,  the  choice  of  the  relaxation  parameter  w  and  the  device  geometry,  i.e.  the  aspect 

ratios  and  drain  separations.  For  example,  a  MOSFET  with  dimensions  W  =  L  =  4  pm,  a  drain 

separation,  d  =  2  pm,  VGS  =  3  V  and  VDS  =  1  V  needs  approximately  1500-2000  SLOR 

-4 

cycles  for  full  convergence.  This  was  for  a  relative  precision,  e  =  10  (see  Fig.  4.6).  Working 
with  double -precision  floating-point  arithmetic  and  a  grid  size  of  50  x  50  nodes  (XMAX  =  50, 
YMAX  =  50),  the  two  dimensional  simulations  of  the  split-drain  magnetic  field  sensor  required 
0.75  to  1  Mbyte  of  memory.  As  for  the  initial  guess,  the  variable  u  was  provided  with  linear 
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values,  i.e. 


u°(M,  N)  = 


UD1  +  uD2 


“  US  /  (YMAX-1) (N-l) 


+  us  (4.36) 


f°r  1  <  N  <  YMAX  and  1  <_  M  <_  XMAX  , 


where  u^  ,  u  ^  and  are  defined  by  (4.19),  (4.20)  and  (4.21)  respectively.  A  value  of  w 
between  1.5  and  1.75  was  found  to  be  an  optimum  range  for  the  nature  of  the  given  initial  guess 
(4.36)  and  for  a  variety  of  drain  separations  and  bias  voltages. 

For  the  Newton -Raphson  iterations  [41],  V  was  provided  with  linear  initial  guess  values, 
i.e. 


V° (M,  N)  = 


V  +  V 
D1  D2 


-  V 


s  /  (YMAX-1) (N-l) 


+  Vc  (4.37) 


for  1  <  N  <  YMAX  and  1  <  M  <  XMAX 


where  V_,  ,  V  „  are  the  voltaees  at  the  drains  and  V  is  the  source  voltaee  (see  Fig.  4.2). 

D1  D2  S  c 

The  CPU  time,  on  the  Amdahl  580,  is  approximately  20-30  seconds  for  the  modeling  of 

the  split -drain  MOSFET's.  At  present,  the  iterations  converge  for  not  too  large  products  of  B 

and  Hall  mobility,  say  tg  ©„  f  0.6. 

Hn 

To  check  if  the  true  solution  was  obtained,  the  computed  results  were  subject  to  the  same 
series  of  tests  described  in  section  3.4.3. 
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Fig.  4.6 


Flow-chart  of  the  solution  procedure 
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4.6  Results  and  Discussion 

The  distributions  of  electric  potential,  current  density  and  surface  charge  are  computed  for 

N-channel  MOSFET's  (uniform  and  split-drain  configurations)  for  a  variety  of  aspect  ratios, 

drain  separations,  operating  voltages  and  field  strengths.  The  magnetic  field  is  perpendicular  to 

the  device  plane.  The  devices  are  in  the  linear  region  of  operation  [42]  with  bulk  doping  density, 
15  -3 

N  .  =  10  cm  ,  oxide  thickness,  t  =  100  nm  and  threshold  voltage,  V_  =  1  V. 

A  ox  l 

Figs.  4.7  and  4.8  compare  equipotential  and  current  lines  in  a  uniform  drain  NM05  [42], 

for  field  strengths  of  1 1  and  3  t.  Here  VDS  =  1  V,  VGg  =  3  V  and  V  g  =  0  V.  In  Fig.  4.9,  we 

see  the  surface  charge  density  for  the  device  operated  close  to  the  pinch -off  region, 

VDS  =  3.2  V,  VGS  =  5  V  and  Vg  =  0  V.  As  the  device  is  operated  closer  to  the  saturation 

region,  a  somewhat  distributed  pinch -off  region  is  found  because  of  the  magnetic  field.  Fig.  4.10 

shows  equipotential  and  current  lines  for  the  device  shown  in  Fig.  4.9.  As  opposed  to  the  results 

obtained  previously  for  the  n-type  silicon  slab  [16,17],  the  potential  lines  tend  to  crowd  at  the 

vicinity  of  the  drain.  Figs.  4.11-4.13  illustrate  equipotential  and  current  lines  for  a  split-drain 

MOSFET  for  different  aspect  ratios,  W/L  and  drain  separations,  d.  The  operating  voltages  are 

V  =  V  =  1  V,  V  =  3  V,  V  =  0  V  and  the  magnetic  field  strength  is  1 1. 

D1S  D2S  GS  S 

The  sensitivity  defined  in  terms  of  the  relative  current  imbalance  in  the  two  drains,  viz. 
6=(I  -I  )/(I  +  I  ),is  studied  for  the  various  MOSFET  geometries  mentioned 

above  [42].  There  is  no  notable  change  in  sensitivity  with  fluctuations  in  the  operating  voltages. 
The  sensitivity  change  becomes  notable  only  when  fluctuations  in  the  gate  bias  are  large,  e.g. 
doubling  the  gate  voltage  leads  to  a  significant  drop  in  sensitivity.  This  is  due  to  the  strong 
dependence  of  the  mobility  of  the  electric  field  in  the  z -direction  (see  Fig.  4.1). 

The  results  suggest  that  the  relative  current  imbalance,  5  ,  depends  more  strongly  on  the 
absolute  values  of  L  or  W  rather  than  the  aspect  ratio  [42,43].  This  is  shown  in  Figs.  4.14  -  4.16 
for  aspect  ratios  W/L  =  25  pm/25  pm,  W/L  =  5  pm/25  pm  and  W/L  =  100  pm/25  pm 
respectively.  It  is  found  that  narrowing  the  channel  down  from  25  pm  to  5  pm  (Fig.  4.15)  leads  to 
only  a  slight  increase  in  sensitivity,  while  making  the  channel  broader  leads  to  a  less  abrupt  carrier 
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deflection  close  to  the  drains,  resulting  in  a  substantial  sensitivity  decrease  (Fig.  4.16). 

From  the  point  of  view  of  keeping  the  device  small,  but  allowing  for  the  minimal  drain 
separation  required  by  the  design  rules,  the  square  geometry  seems  to  be  favourable  [43].  A  circuit 
configuration  matched  to  the  square  device  (W  =  L  =25  pm,  d  -  5pm)  is  shown  in  Fig.  4.17.  To 
achieve  the  maximum  attainable  sensitivities  the  drain  voltages  V  ,  V  are  put  at  the  same 

L)  -L  o  JJZ  o 

potential  in  order  to  avoid  any  shunting  effects  between  the  drains  and  the  device  is  operated  at 

low  gate  voltages,  e.g.  V  <  V  <  5  V. 
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Fig.  4.7  Equipotential  and  current  lines  for 
W  =  L  =  4  }im,  B  =  1  t 


Fig.  4.8  Equipotential  and  current  lines  for 


W  -  L  =  4  jim,  b 


3  t 
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Surface  charge  density  for  W  =  L  =  4  g  =  1  t 


Fig.  4.9 
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Fig.  4.10 


Equipotential  and  current  lines  for  W  =  L  =  4  pm, 
B  =  1  t 
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Fig.  4.11 


Equipotential  and  current  lines  for  W/L  =  4  pm/4  ^ 
d  =  1.4  pm,  Brit 
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Fig.  4.12  Equipotential  and  current  lines  for 

W/L  =  22  pm/14  pm,  d  =  7.63  pm,  B  =  1  t 


I 


Fig.  4.13 


Equipotential  and  current  lines  for 

W/L  =  14  pm/22  pm,  d  =  4.9  pm,  B  =  1  t 


69 


Fig.  4.14  Relative  current  imbalance,  6  as  a  function  of  magnetic 
field,  B  for  a  square  channel  NMOS,  W  =  L  =  25  ym, 
for  different  relative  drain  separations  : 

d/W  =  0.05  (curve  a),  d/W  =  0.2  (curve  b) ,  d/W  =  0.6  (curve  c) 


Fig.  4.15  Same  as  Fig.  4.14  but  W/L  =  5  ym/25  Pm 


«[*]  A 

15 


10 


-10 


-15 


Fig.  4.16  Same  as  Fig.  4.14  but  W/L  =  100  ym/25  urn 


Fig.  4.17  Schematic  of  the  square  channel  device  in  circuit 


5.  CONCLUSION  AND  OUTLOOK 

Semiconductor  devices  subject  to  an  external  magnetic  field  have  been  modeled 
numerically  using  standard  methods  [3-5,28,44]  but  incorporating  the  magnetic  field  into  the 
current  density  equations  and  the  boundary  conditions.  To  obtain  a  stable  solution  method,  the 
generalisation  of  the  well-known  Scharfetter-Gummel  scheme  has  been  used  in  the  case  of  two 
dimensions  and  nonzero  magnetic  field.  The  results  show  the  important  point  that  the  two 
dimensional  distributions  of  potential  and  current  density  cannot  be  easily  described  in  terms  of 
the  simple  analytical  models  such  as  Hall  voltage,  Lorentz  deflection  or  magnetoconcentration  ; 
this  is  due  to  the  intricate  combination  of  the  magnetic  field  and  the  boundary  conditions.  Only  in 
certain  limits  of  the  device  geometry  can  the  action  of  the  field  be  visualised  in  terms  of  these 
simple  analytical  models. 

Presently,  the  algorithm  is  being  extended  to  p+-i-n+  devices  and  triple-drain 
MOSFET's.  In  the  p  +  -i-n+  device,  split -contacts  are  incorporated  and  the  sensitivity  of  the 
device  to  the  magnetic  field  is  measured  as  the  relative  current  imbalance  between  the  contacts 
[36].  It  has  been  found  that,  if  the  diffusion  lengths  of  the  carriers  are  adjusted  suitably  to  the 
width  and  length  of  the  p+  -i-n+  device,  sensitivity  figures  superior  to  those  of  Hall  plate  devices 
are  achieved.  With  split -drain  MOSFET's,  by  introducing  an  additional  drain  in  the  'split'  region, 
the  sensitivity  figures  obtained  are  substantially  higher  than  those  of  dual-drain  devices  [45],  The 
sensitivity  of  the  triple -drain  MOSFET's  is  taken  to  be  the  relative  current  imbalance  between  the 
outer  two  drains. 

With  these  results,  future  work  may  extend  the  numerical  algorithm  to  magnetotransistors 
fora  variety  of  device  configurations  and  investigating  the  signal -to -noise  ratio  limitations  on  the 
sensitivity  of  the  devices  modeled  so  far.  In  low  frequency  applications  (e.g.  60  Hz),  the  dominant 
type  of  noise  in  MOSFET’s  is  flicker  (1/f)  noise.  Efforts  should  be  concentrated  on  being  able  to 
include  the  magnetic -field  effects  in  one  of  the  existing  circuit  simulation  programs,  e.g.  SPICE, 
thus  allowing  the  evaluation  and  testing  of  the  different  sensor  models  for  future  computer  aided 
design  of  monolithic  integrated  magnetic -field  sensors  that  could  include  digital  signal  processing 
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such  as  sampling  and  A/D  conversion.  This  would  also  allow  the  optimisation  of  the  design  rules 
such  that  'stray'  magnetic  fields  have  little  or  no  influence  on  device  characteristics  when  operated 
in  a  magnetic  environment. 

Integrated  Si  MFS  share  problems  of  noise,  temperature -coefficient,  offset,  nonlinearity 
and  frequency -dependence  with  comparable  nonmagnetic  devices  and  analog  IC's  of  the  same 
technology,  e.g.  the  offset  due  to  nonideal  device  geometry  and  the  piezoresistive  effect  of  Si. 
Some  of  those  problems  can  possibly  be  overcome  by  using  procedures  known  from  analog  IC 
design  and  careful  process  technology  [1]. 
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APPENDIX  A 


Linear  Interpolation  of  the  Variables  p,  n  and  \p. 


The  interpolation  of  the  variables  p,  n,  ip  and  the  simplification 
of  the  Bernoulli  coefficients  used  in  obtaining  (3.30)  is  explained 
in  this  appendix.  With  reference  to  Fig.  3.3, 
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^BA  =  *B  _  ♦a  =  1/2  (*3  '  V  ' 


^CB  =  “  ’"b  =  1/2  (*4  "  *?>  - 


4  r2 
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The  Bernoulli  coefficients  in  the  expanded  version  of 
eqn.  (3.30)  can  be  simplified  as  follows; 
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Similarly, 
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At  this  stage  eqn.  (3.30)  now  reads 
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By  linearly  interpolating  nA,  nB,  nc  and  nD  in  the  same 
fashion  ip^,  and  was  interpolated,  we  obtain 
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Hence  (A.l)  can  be  further  simplified  to, 
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In  obtaining  (A. 2)  the  following  simplification  was  used 
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The  hole  continuity  equation  was  simplified  on  the 
same  basis  of  interpolation. 


(A. 2) 


A  better  interpolation  technique  for  the  variables  p,  n  and  \p 
maintaining  consistency  with  Scharfetter-Gummel  scheme  is 
described  below. 
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♦a  =  I  (h  +  *2)  ' 

*B  =  |  (V>i  +  V)  ' 

*0  =  1  (*3  +  <P4)  , 

I'D  =  I  (*i  +  V*  • 


Equating  the  integrals  obtained 
over  s  (first  between  i  and  A1 , 
after  simplifications. 


when  eqn.  (3.21)  was  integrated 
then  between  i  and  A)  we  obtain 


In  a  similar  fashion  n  ,  nc  and  nD  are  obtained.  The  same 
treatment  holds  for  the  simplification  of  the  hole  continuity 
equation. 


APPENDIX  B 


Matrix  and  Vector  Elements  for  the  Slab  Boundaries 


The  matrix  and  vector  elements  for  the  floating  boundaries  of  the 
slab  (see  Fig.  3.2),  given  by  eqns.  (3.44)  and  (3.46),  are  described 
in  this  appendix.  TO  avoid  the  lengthy  expressions  contained  in  the 
matrix  and  vector  elements,  we  define  the  following  ; 
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Here,  p° (M,  N) ,  n° (M,  N)  and  ^°(M,  N)  are  initial  guess 
values  specified  at  every  meshpoint. 


The  matrix  elements  are: 
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The  vector  elements  are: 

kl  =  A1  *0  +  A0  Tp  V  .  k2  =  *0  +  A0  Tn  Y  , 
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For  the  floating  boundary  (see  Fig.  3 . 2) , 
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For  the  floating  boundary  (^T^)  (see  Fig.  3.2)  , 
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APPENDIX  C 


Matrix  and  Vector  Elements  for  the  Bulk 


In  this  appendix,  the  matrix  and  vector  elements  related  to  the 
coefficients  of  the  continuity  and  Poisson's  equations  given  by 
eqn.  (3.50),  are  described.  To  avoid  the  lengthy  expressions 
contained  in  the  elements,  the  following  definitions  are  made  : 
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ABSTRACT 

We  present  the  first  two-dimensional  numerical 
solutions  of  the  coupled  partial  differential  equa¬ 
tions  governing  carrier  drift,  diffusion,  genera¬ 
tion,  recombination  and  electric  potential  in  the 
bulk  of  a  semiconductor  material  in  the  presence  of 
a  magnetic  field.  Our  solutions  were  obtained  for 
a  uniformly  doped,  quadratic  semiconductor  slab 
subject  to  boundary  conditions  which  arise  from 
two  ohmic  contacts  supplying  the  device  current. 


I  -+  -*• 

-  7  •  J  -  R  =  0, 

q  n 


(7)^  *  - -3.  (p  -  n  +  N) , 


(3) 

(4) 


(5) 


INTRODUCTION 

In  the  presence  of  a  magnetic  field,  numerical 
modeling  of  a  semiconductor  device  is  required 
because  of  the  complexity  of  the  situation;  in  the 
bulk,  the  magnetic  field,  electric  field,  drift 
current,  diffusion,  generation  and  recombination  of 
the  carriers  are  all  coupled  to  one  another.  More¬ 
over,  the  actual  bulk  spatial  distributions  of 
these  quantities  in  the  device  are  determined  by 
the  boundary  conditions  specified  due  to  the 
various  (open  or  ohmic)  device  contacts. 

BASIC  EQUATIONS 

In  this  paper,  we  present  two-dimensional 
numerical  solutions  of  the  following  system  of 
coupled  partial  differential  equations  [1,2]: 


(pn  -  n;.2) 

E.-E  E  -E. 

t  [p  +  n  exp (~ j/T  ")]  +  t  [n  +  n  exp(-[  X)]. 
n  l  XI  p  i  kl 


Here,  Jn  and  Jp  represent  electron  and  hole  current 
densities,  un  and  Up  the  drift  mobilities,  and 
Up  the  Hall  mobilities,  n  and  p  the  electron  and 
hole  concentrations,  D  and  Dp  the  diffusion 
constants,  E  the  electric  field,  B  the  magnetic 
induction,  ^  the  electrostatic  potential,  7  the 
nabla  operator,  N  the  net  doping  density,  q  denotes 
the  elementary  charge,  R  is  the  net  bulk  recombina¬ 
tion  rate  (SRH  type),  n^  the  intrinsic  concentra¬ 
tion,  xn  and  ip  denote  the  bulk  carrier  lifetimes, 
Et  the  energy  level  of  trap  centres,  E^  the 
intrinsic  Fermi  level,  k  the  Boltzmann’s  constant 
and  T  the  absolute  temperature. 
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Equations  3  and  4  represent  the  continuity 
equations  for  the  stationary  case  and  equation  5 
represents  Poisson's  equation. 
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In  this  paper  we  consider  a  quadratic  semi¬ 
conductor  slab  with  open  boundaries  along  two 
opposite  faces  and  ohmic  contacts  along  the  other 
two  opposite  faces  of  the  slab.  The  magnetic  field 
is  assumed  to  be  perpendicular  to  the  device  sur¬ 
face:  B  *  (0,  0,  Bz),  so  that  all  distributions 
in  the  device  depend  only  on  x  and  y. 

For  the  mobilities  we  used  [3] 
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Fig.  1.  Equipotential  curves  for  the  n-type 
silicon  device,  -10lfccm  V  =0.1  volt, 
L  -  W  =  lOOu  and  B  =  1  Tesla.  a 
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where  Nj  denotes  the  total  impurity  concentration. 
For  the  Hall  mobility,  we  used  the  relation,  u*  = 
ru  where  we  adopted  for  r  the  value  1.92  [4], 


BOUNDARY  CONDITIONS 

Along  the  ohmic  contacts,  we  assumed  [2,3] 

2 

pn  -  ni  .  (9) 

Furthermore  along  the  ohmic  contacts  we  obtain  for 
the  electrostatic  potential,  from  the  condition  of 
charge  neutrality  [2,3] 

*  -  va  +  “  Sinh"1(2^7)’  (10) 


where  Va  denotes  the  applied  voltage. 

Along  the  open  boundaries,  we  have  used  the 
conditions  [2,3]: 
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In  computing  the  equations  (11), 
equations  (1)  and  (2),  which  contains 
field. 


we  used  the 
the  magnetic 


Because  of  the  non-zero  magnetic  field  along 
the  open  boundary,  we  have  chosen  the  boundary 
condition 


~  (p  -  n  +  N)  -  0. 
3  x 


(12) 


This  condition  is  considered  somewhat  less  restric¬ 
tive  than  proposed  by  Kurata  [2]. 


NUMERICAL  PROCEDURE 

Following  the  technique  given  in  References 
[2,3,5]  we  have  put  the  equations  (1)  -  (12)  above 
in  finite  difference  form.  For  the  continuity 
equations,  (3)  and  (4),  which  now  include  the 
magnetic  field,  we  used  a  discretization  technique 
similar  to  a  generalization  of  the  Scharfetter- 
Gummel  scheme  discussed  recently  by  Rudan  and 
Guerrieri  [6]. 

Our  iteration  procedure  is  based  on  (i)  the 
linearization  of  the  equations  by  applying 
Newton's  iteration  principle,  and  (ii)  the 
simultaneous  solution  of  Poisson’s  and  the  con¬ 
tinuity  equations  ("full  coupling")  [2], 


ig.  2. 


Current  lines  corresponding  to  Fig.  1. 
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Fig.  3.  Equipotential  curves  for  the  p-type 
silicon  device,  N___=  lO^cm  Vg  =  0.1  volt, 
L  =  W  =  100_  and  I  =  10  Tesla. 


RESULTS 

We  computed  the  electric  potential  and  the 
current  density  in  a  square  silicon  slab  with  ohmic 
contacts  for  magnetic  fields  up  to  10  Tesla,  and 
with  a  doping  level  of  lO^cnT^  for  both  p  and  n 
type  material. 

Figures  1—4  show  plots  of  the  equipotential 
curves  and  the  current  lines.  Under  the  conditions 
that  apply  to  our  computations,  magnetoconcentra¬ 
tion  is  found  to  be  insignificant.  Consequently, 
we  did  not  make  plots  of  the  carrier  densities. 

Also,  our  results  depend  very  little  on  the  linear 
dimensions  of  the  square  slab. 

Due  to  the  boundary  conditions,  equipotential 
lines  close  to  the  metal  have  to  be  parallel  to  the 
boundary.  Thus  carriers  are  deflected  in  this  area. 
Close  to  the  open  boundary,  however,  current  lines 
are  parallel  to  the  boundary;  a  Hall  electric  field 
is  built  up  in  this  region. 

Far  from  the  boundaries,  however,  the  opera¬ 
tion  of  the  device  can  be  visualized  neither  in 
terms  of  Lorentz  deflection  nor  Hall  voltage;  both 
of  these  mechanisms  are  involved  in  a  complex  way. 

The  deviation  of  our  results  from  the  limit 
W/L  -*■  0  (Vinal's  case  [7])  is  demonstrated  by  the 
variation  of  the  electric  field  component  ratio 
Ex / Ey  (see  Figs.  5-6)  across  the  device.  In  the 
limit  W/L  -*■  0  this  ratio  is  a  constant. 


DISCUSSION 

To  the  best  of  our  knowledge,  we  are  presen¬ 
ting  the  first  two-dimensional  modeling  of  a 
magnetic-field-sensitive  semiconductor  device.  Our 
results  show  the  important  point  that  the  two- 
dimensional  distributions  of  current  density  and 
potential  cannot  be  easily  described  in  terms  of 
simple  models  such  as  Hall  voltage  or  Lorentz 
deflection;  this  is  due  to  the  intricate  combina¬ 
tion  of  the  magnetic  field  and  the  boundary 
conditions. 

Thus  our  results  lead  to  an  accurate  predic¬ 
tion  of  the  Hall  voltage  even  under  extreme  condi¬ 
tions  in  terms  of  the  geometries  and  the  magnetic 
fields  such  that  the  Hall  voltage  is  a  nonlinear 
function  of  the  magnetic  field  strength. 

Further  work  concerning  magnetic  field 
sensors  that  comprise  intrinsic  semiconductor 
materials  and  p-n  junctions  is  in  progress. 
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Fig.  4.  Current  lines  correpsonding  to  Fig.  3. 
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Fig.  5.  Electric  field  ratio  Ex/Ey  across  the 
middle  (y  =  L/2)  of  device  corresponding  to  Fig.  1. 
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Fig.  6.  Electric  field  ratio  Ex/Ey  across  the 
middle  (y  =  L/2)  of  device  corresponding  to  Fig.  3. 
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Two-Dimensional  Numerical  Analysis  of  a 
Silicon  Magnetic  Field  Sensor 
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Abstract -W e  present  two-dimensional  numerical  solutions  of  the 
coupled,  nonlinear,  partial  differential  equations  governing  the  electric 
potential,  carrier  drift,  diffusion,  generation,  and  recombination  in  a 
finite  semiconductor  slab  in  the  presence  of  a  magnetic  field.  This 
enables  device  modeling  for  general  geometries,  doping  levels,  and 
injection  conditions,  where  the  effect  of  the  magnetic  field  cannot 
be  expressed  simply  in  terms  of  Hall  voltage,  Lorentz  deflection,  or 
magnetoconcentration. 


Introduction 

Currently,  we  see  rapid  growth  in  the  area  of  magnetic-field- 
sensitive  silicon  microtransducers  (see,  e.g.,  [  1  ]  — [4]  and  refer¬ 
ences  therein).  The  theoretical  understanding  of  the  electronic 
mechanisms  underlying  these  devices,  let  alone  device  model¬ 
ing,  has  not  kept  up  with  this  rapid  growth.  This  situation 
eventually  has  led  to  some  dispute  [5] ,  [6] . 

Device  modeling  in  the  absence  of  a  magnetic  field  is  a 
highly  developed  area  of  research.  Many  devices  can  only  be 
modeled  by  numerical  solution  of  the  well-known  nonlinear 
system  of  partial  differential  equations  comprising  the  conti¬ 
nuity  equations  and  Poisson’s  equation  [7]  —[10]. 

In  the  presence  of  a  magnetic  field,  the  numerical  modeling 
of  semiconductor  devices  is  even  more  required  for  the  follow¬ 
ing  reasons: 

i)  In  the  bulk,  the  magnetic  field  causes  a  sequence  of  ef¬ 
fects  to  take  place:  changes  of  the  carrier  motion,  the  carrier 
distributions,  and  the  electrostatic  potential.  Only  in  certain 
limits  the  action  of  the  field  can  be  visualized  simply  in  terms 
of,  e.g.,  Hall  voltage  or  carrier  deflection. 

ii)  Any  realistic  device  modeling  has  to  account  for  the 
boundary  conditions,  such  as  “ohmic”  or  “oxide”  contacts 
[10],  The  magnetic  field  introduces  some  asymmetry  into  the 
current  equations,  while  the  boundary  geometry  retains  any 
zero-field  symmetry;  this  makes  the  problem  even  more 
difficult. 
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Basic  Equations 

Because  of  the  Lorentz  force,  the  current  equations  in  a 
semiconductor,  in  the  presence  of  a  magnetic  field  read  as  fol¬ 
lows  [11]: 

-*  1  .*  -* 

Jn  =  1 +  tf>nVn) 

1  +  (MnS)2 

+  P%B  X  (qn„nE  +  qDnVn)]  (1) 

-»  1  4  -• 

jp  ~ - 77T7  [(qHppE  -  qDpS/p) 

P  1  +(MpB)2  P  P 

-  PpB  X  ( qPppE  -  qDp\Jp )] .  (2) 

Here,  B  is  the  magnetic  induction,  Jn  and  Jp  denote  the  elec¬ 
tron  and  hole  current  densities,  pn  and  pp  the  drift  mobilities, 
p*  and  Up  the  Hall  mobilities,  n  and  p  the  electron  and  hole 
concentrations,  Dn  and  Dp  the  diffusion  coefficients;  q  de¬ 
notes  the  elementary  charge,  E  the  electric  field,  and  V  the 
nabla  operator.  The  continuity  equations  (stationary  case) 
and  Poisson’s  equation  are  as  usual,  viz., 


1  -»  ■* 

-  V  •  /„  -  R  =  0 
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l  -►  -* 
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q  p 
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By  R  we  denote  the  net  bulk  recombination  rate,  N  the  net 
doping  density,  e  the  dielectric  constant,  and  \ p  the  electro¬ 
static  potential.  For  R,  we  use  the  Schockley-Read-Hall 
model,  viz., 

_ ( pn  ~  nj ) _ 

t„  [p  +  rii  exp  (  ^f-)]  +  TP  ["  +  ni  exP  (  \T  ' 

(6) 

Here  r„  and  Tp  are  the  bulk  carrier  lifetimes,  Et  is  the  energy 
level  of  trap  centers,  the  intrinsic  Fermi  level,  k  is  the  Boltz¬ 
mann’s  constant,  T  the  absolute  temperature,  and  nt  denotes 
the  intrinsic  concentration.  Following  common  practice  [7], 
in  our  computations  we  assumed  Et  =  Ej. 

Equations  (1)  to  (6)  constitute  the  system  to  be  solved  sub¬ 
ject  to  appropriate  boundary  conditions.  In  this  brief,  we  dis¬ 
cuss  only  a  rectangular  semiconductor  slab  with  the  magnetic 
field  perpendicular  to  the  device  surface,  B  -  (0,  0,  Bz),  as 
shown  in  Fig.  1 ;  then  all  distributions  in  the  device  depend 
only  on  x  and  y. 

For  the  mobilities  we  used  [  12] 
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where  Nf  denotes  the  total  impurity  concentration.  For  the 
Hall  mobility,  we  used  the  relation,  p*  =  rp  where  we  adopted 
for  r  the  value  1 .92  [  1 3 1 . 
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Fig.  1.  Basic  device  geometry:  semiconductor  slab  of  length  L  and 
width  W,  with  a  magnetic  field  6  perpendicular  to  the  drawing  plane. 
The  applied  voltage  is  Va. 


Boundary  Conditions 

We  take  ohmic  contacts  on  two  opposite  faces  of  the  device, 
as  shown  in  Fig.  1.  Along  these  contact  boundaries,  we  as¬ 
sume  [7] , [9] 

pn  =  nf .  (9) 

Moreover,  charge  neutrality  leads  to  the  condition  [7] ,  [9] 

\p  =  Va  +  —  sinh'1  ( - \  (10) 

for  the  potential  \J/.  Here,  Va  denotes  the  applied  voltage  and 
kT/q  equals  25.7  mV  at  room  temperature. 

Along  the  other,  “floating  boundaries,”  we  maintained  the 
current  conditions  [7] ,  (9] 

Jnx-0  Ola) 

Jpx  -  0.  Olb) 

Computing  the  boundary  conditions  (1 1),  we  used  (1)  and  (2) 
which  contain  the  magnetic  field. 

As  for  the  space  charge  condition  along  the  floating  bound¬ 
ary,  we  have  adopted  the  condition 


Fig.  2.  Equipotentia!  lines  and  current  lines  for  a  square  device,  L  = 
W  =  100  Mm,  n-type  silicon,  doping  level  10 16  cm-3,  Bz  =  1  T. 
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Fig.  3.  (a)  Same  as  Fig.  2,  but  L  =  4 IV,  Bz  =  1  T.  (b)  Same  as  Fig.  2, 

but  L  =  4W,  Bz  =  2  T. 


— (p  -  n  +  N)  =  0.  (12) 

bx2 

Equation  (12)  is  less  restrictive  than  the  condition  of  zero 
space  charge  variation  proposed  by  Kurata  [9] .  The  choice  of 
(12)  was  required  because  of  the  asymmetry  introduced  to  our 
problem  by  the  magnetic  field. 


Numerical  Procedure 

We  have  set  up  ( 1  >— ( 12)  in  finite  difference  form  following 
the  technique  of  [9]  and  [12].  Moreover,  for  the  continuity 
equations  (3)  and  (4),  which  include  the  field  B  through  (1) 
and  (2),  we  have  used  a  discretization  technique  similar  to  the 
generalization  of  the  Scharfetter-Gummel  scheme  discussed 
recently  by  Rudan  and  Guerrieri  [  14] . 


I 


Our  iteration  procedure  is  based  on  i)  the  linearization  of  the 
equations  by  applying  Newton’s  iteration  principle,  and  ii)  the 
simultaneous  solution  of  Poisson’s  and  the  continuity  equa¬ 
tions  (“full  coupling”)  [9] . 

Results 

To  demonstrate  the  dependence  of  the  solutions  of  the  de¬ 
vice  geometry,  we  computed  the  electric  potential  and  the  cur¬ 
rent  density  for  a  variety  of  ratios  W/L  of  width  to  length  of 
the  semiconductor  slab.  As  an  example,  we  chose  n-type  sili¬ 
con  with  a  doping  level  1016  cm"3;  the  magnetic  field  varied 
from  0  to  2  T;  the  applied  voltage  Va  was  0.1  V. 

Plots  of  the  equipotential  lines  and  current  lines  are  shown 
for  W/L  =  1  (Fig.  2),  W/L  =  ±  (Fig.  3),  and  W/L  =  4  (Fig.  4(a)). 
The  results  strongly  depend  on  the  ratio  W/L,  but  depend  very 
little  on  the  absolute  values  of  W  or  L,  provided  these  linear 
dimensions  are  not  too  small  (larger  than  say,  1  /Ltm). 

In  the  long-slab  limit  {W  «L,  see  Fig.  3),  our  results  resem¬ 
ble  those  of  Vinal’s  magnetic-field-sensor  model  [15],  [16], 
i.e.,  there  is  a  constant  Hall  field  Ex,  and  the  current  lines  are 
parallel  to  the  slab  boundaries,  except  for  regions  close  to  the 
ohmic  contacts.  In  the  opposite  limit  of  a  short  device  with 
broad  contacts  (1 V»L,  see  Fig.  4(a))  carrier  deflection  pre¬ 
vails  in  most  of  the  device  area,  viz.,  there  are  horizontal  equi¬ 
potential  lines  and  skew  current  lines  showing  the  usual  Hall 
angle,  as  discussed  by  Zieren  [17] . 

In  the  general  case  between  these  extremes  (e.g.,  the  square- 
shaped  device,  see  Fig.  2),  the  operation  of  the  device  can  be 
visualized  neither  in  terms  of  Lorentz  deflection  nor  Hall  volt¬ 
age;  both  these  mechanisms  are  involved  in  a  complex  way.  It 
turns  out  that  magnetoconcentration  is  insignificant  under  the 
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Fig.  4.  (a)  Same  as  Fig.  2,  but  L  =  W/4,  Bz  *  1  T.  fb)  Equipotential 
lines  and  current  lines  for  a  p-type  silicon  device,  doping  level  10 16 
cm"3,  L  -  W/4,  Bz-  2  T. 


circumstances  considered  in  this  paper.  However,  we  found 
magnetoconcentration  to  be  significant  in  other  cases,  e.g., 
under  intrinsic  conditions. 


The  deviation  from,  say,  Vinal’s  limit  W/L  -*•  0,  can  be  dem¬ 
onstrated  in  terms  of  the  electric  field  component  ratio  Ex/Ey 
or  the  current  deflection  Jx/Jy  across  the  device.  We  note  that 
both  ratios  Ex/Ey  and  Jx/Jy  have  the  maximum  value  of  tan 
Qh  -  H*BZ  where  0#  is  the  Hall  angle  (see  Figs.  5  and  6). 

We  have  furthermore  computed  the  case  of  p-type  silicon 
with  1016  cm"3  doping  level  for  fields  up  to  5  T.  Analogous 
results  are  obtained;  the  magnetic  field  effects  are  somewhat 
weaker  because  of  the  lower  Hall  mobility  for  holes  (see 
Fig.  4(b)). 


Discussion 

To  the  best  of  our  knowledge,  in  this  brief  we  are  reporting 
the  first  numerical  modeling  of  a  magnetic-field-sensitive  semi¬ 
conductor  device.  We  have  used  standard  methods  [7] -[10], 
[12],  [14]  appropriately  adapted  to  the  problem  with  exte¬ 
rior  magnetic  field.  The  crucial  point  in  modeling  such  a  de¬ 
vice  is  the  intricate  combination  of  the  magnetic  field  and  the 
boundary  conditions;  this  may  result  in  the  failure  of  simpler 
analytical  models  based  solely  on,  e.g.,  the  Hall  effect  or  the 
Lorentz  deflection  to  describe  these  devices. 

Our  present  iteration  converges  for  not  too  large  products 
of  Bz  and  Hall  mobility,  say  tan  0//  <  0.3.  For  all  numerical 
computations  reported  in  this  brief,  the  use  of  a  uniform  grid 
was  found  to  be  sufficient.  Work  is  in  progress  concerning  de¬ 
vices  in  the  presence  of  an  exterior  magnetic  field  under  intrin¬ 
sic  conditions. 
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Fig.  6.  Current  deflection  JxUy  across  the  middle  ( y  =  LI 2)  of  the  de¬ 
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APPENDIX  F 

NUMERICAL  MODELING  OF  SILICON  MAGNETIC  FIELD  SENSORS: 
MAGNETOCONCENTRATION  EFFECTS  IN  SPLIT-METAL-CONTACT  DEVICES 
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Abstract :  The  two-dimensional  distributions  of  the 
electric  potential,  the  electron  concentration  and  the 
hole  concentration  in  a  silicon  slab  exposed  to  a  mag¬ 
netic  field  have  been  computed  numerically.  Generalizing 
the  well-known  Scharfetter-Guramel  scheme  to  the  case 
of  two  dimensions  and  nonzero  magnetic  field,  we  have 
employed  a  finite-difference  technique.  In  case  of 
Hall  plates,  where  space  charge  is  negligible,  our 
results  are  in  support  of  previous  results  obtained 
by  analytical  models  or  by  conformal-mapping  technique. 
In  intrinsic  or  closely  intrinsic  silicon,  our  results 
show  both  magnetoconcentration  and  space-charge  effects; 
these  can  be  exploited  in  the  design  of  split-metal- 
contact  magnetic  field  sensors. 

Introduction 


Here  n  and  p  are  the  electron^and  ho^e  concentrations, 

Y  is  the  electric  potential,  Jn  and  Jp  are  the  electron 
and  hole  current  densities,  q  is  the  elementary  charge, 
c  is  the  dielectric  constant,  R  is  the  net  bulk  recom¬ 
bination  rate  and  N  the  ionized  net  doping. 


When  the  magnetic  field  is  perpendicular  to 
the  current  densities,  the  transport  equations  read  [4], 


1  -+ 

- TTT  Kqu  n  E  +  qD  vn) 

1  +  (U  B)2  n  n 
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+  u  B  X  (qy  nF  +  qD  7n) ] , 
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A  large  number  of  new  magnetic-field-sensitive 
silicon  devices  have  been  realized  recently  (see 
[1]  and  references  therein).  The  theoretical  under¬ 
standing  of  the  electronic  mechanisms  underlying 
these  devices,  let  alone  device  modeling,  has  not 
kept  pace  with  this  rapid  growth.  This  situation 
has  eventually  led  to  some  dispute  [1], 

Numerical  semiconductor  device  modeling,  a  well- 
established  art  for  zero  magnetic  field,  is  even 
more  desirable  for  magnetic  sensors: 

The  basic  magnetic  effects  (Hall  voltage, 

Lorentz  deflection,  magnetoconcentration)  are 
two-dimensional  and  cannot  be  simply  super¬ 
imposed  to  describe  the  device.  Moreover,  the 
role  of  the  boundary  conditions  is  crucial, 
since  the  magnetic-field  effects  are  particularly 
strong  close  to  the  device  boundaries. 

The  problem  in  question  is  made  more  difficult 
by  the  fact  that  the  magnetic  field  breaks  some  of 
the  symmetries  exploited  in  the  modeling  of  zero- 
field  semiconductor  devices.  Apart  from  some  early 
work  in  terms  of  the  Hall  plate  model  [2],  numerical 
results  on  silicon  magnetic  field  sensors  have 
become  available  only  recently  [1J.  While  our 
previous  note  emphasizes  geometrical  aspects  of  semi¬ 
conductor  slabs  with  metal  contacts  in  a  magnetic 
field,  this  paper  presents  new  results,  focused  on 
the  magnetoconcentration  effect  and  split-contact 
(current  imbalance)  applications.  For  a  full  presen¬ 
tation  of  the  underlying  numerical  methods  we  refer 
to  [3]. 


Equations  and  Boundary  Conditions 

In  a  stationary  magnetic  field,  the  continuity 
equations  for  electrons  and  holes  and  the  Poisson 
equation  are 


■  — V-  2  I(qv^  -  <vp> 

P  1  +  (y  B)  P  P 
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A  *+■  ^ 

-  u  B  x  (qu  pE  -  qD  Vp) ] . 

p  ^  pr  p 
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Here,  B  is  the  magnetic  induction,  yn  and  yp  are  the 

drift  mobilities,  u$  and  y£  the  Hall  mobilities, 

D  and  D  the  diffusion  coefficients, 
n  p 


For  the  bulk  recombination  rate  we  use  the 
Schockley-Read-Hall  model  [5],  [6].  For  the  drift 
mobilities  we  use  well-known  empirical  expressions 
[7].  For  the  Hall  mobilities  we  use  the  approxima¬ 
tion  [8]  y*  *  ry,  with  r  given  the  value  1.92  [9]. 

In  this  paper,  we  consider  semiconductor  slabs 
of  rectangular  geometry  with  a  magnetic  field  perpen¬ 
dicular  to  the  device  surface:  B  *  (0,  0,  B  ) 

(see  Fig.  1).  Then  the  distributions  of  carriers 
and  potential  in  the  device  depend  only  on  x  and  y; 
we  solve  a  two-dimensional  problem. 
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Fig.  1  Device  geometry: 
semiconductor  slalj  with 
a  magnetic  field  B  per¬ 
pendicular  to  the  draw¬ 
ing  plane.  The  applied 

voltage  is  V  . 

a 


Fig.  2  Geometry  of  a 
silicon  magnetic  field 
sensor  with  a  split-metal 
contact.  W  and  L  denote 
the  device  width  and 
length  respectively. 


We  assume  metal  contacts  on  two  opposite  faces 
of  the  device  as  shown  in  Fig.  1.  Along  these 
"fixed"  boundaries,  the  values  of  n,  p,  and  ¥  are 
given  by  infinite  recombination  rate  (6),  charge 
neutrality  (7),  and  the  applied  voltage  V  (8) 

3 


'h>-  x 
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[5],  [6],  Using  (kT/q)*  Dn/pn,  we  write 


P  •  n  -  n 

(6) 

p  -  n  +  N  * 

0 

(7) 

ip  *  V  +  — 
a  q 

slnh  (^r>. 

(8) 

conditions  on  the  equipotential  and  current  lines 
can  be  described  as  follows: 

(i)  Close  to  the  metal  boundaries  the  equi¬ 
potential  lines  are  parallel  to  the 
boundaries . 

(ii)  Close  to  the  open  boundaries  the  current 
lines  are  parallel  to  the  boundaries. 


Along  the  other,  insulated  "floating"  boundaries 
we  maintain  the  usual  zero-normal-current-density 
conditions  (see  Fig.  1) 


J 

nx 
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J 

px 
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For  the  third  condition  required  along  the  floating 
boundaries,  we  use 


2-z  (p  -  n  +  N)  -  0.  (11) 

3x 


While  this  condition  is  weaker  than  the  conditions 
commonly  used  at  floating  boundaries,  it  still 
ensures  that  the  space  charge  is  a  smooth  function 
of  the  spatial  coordinates. 

Numerical  Procedure 


The  partial  differential  equations  (1)  -  (5) 
and  the  boundary  conditions  (6)  -  (11)  were  dis¬ 
cretized  for  obtaining  numerical  solutions.  In 
view  of  the  rectangular  geometry  of  our  model 
devices  we  used  for  discretization  the  finite- 
difference  technique  [5]  implemented  on  a  non- 
uniform,  rectangular  grid.  This  technique,  which 
allowed  us  to  use  the  SLOR  method  [5],  [10]  achieves 
reasonably  fast  convergence  even  when  the  equations 

(I)  -  (5)  were  solved  simultanously  by  a  "coupled" 

[II]  or  "simultaneous"  [5]  algorithm. 

After  discretization,  equations  (1)  -  (5)  were 
linearized  by  a  Newton-iteration  scheme  [5].  The 
linear  systems  thus  obtained  were,  however,  in  many 
cases  of  practical  interest  unsuited  for  an  intera- 
tive  solution  because  of  the  numerical  values  the 
out-of-diagonal  elements  took  on.  We  solved  this 
problem  [3]  by  generalizing  the  well-known  Scharfetter- 
Gummel  scheme  [12]  for  the  case  of  two-dimensions  [11], 
[13]  and  nonzero  magnetic  field. 

Using  a  grid  with  50  x  50  nodes  and  working 
with  double-precision  arithmetic  the  computation 
of  one  solution  required  approximately  2  Mbytes  of 
memory  and  50  -  200  seconds  of  CPU  time  on  our 
Amdahl  580  computer.  For  every  Newton  cycle,  15  - 
30  SLOR  iterations  were  implemented.  For  a  relative 
precision  of  10“^  typically  8-20  Newton  cycles 
were  required.  In  the  SLOR  iteration,  we  employed 
a  relaxation  factor  between  0.1  and  1.9  depending 
on  the  device  parameters  and  the  initial  guess 
value . 


Results 


Typical  results  obtained  in  silicon  of  various 
doping  levels  for  the  electric  potential,  the  current 
density,  the  electron  and  hole  concentrations  are 
shown  in  Figs.  3-5.  The  effect  of  the  device 
geometry  has  been  presented  elsewhere  [1]. 

Fig.  3  and  Fig.  5a  show  plots  of  equipotential 
and  current  lines.  The  effects  of  the  boundary 


In  the  case  of  doped  material,  the  results 
(see  Fig.  3)  agree  well  with  those  obtained  earlier 
by  the  Hall-plate  model  [2],  where  the  Poisson 
equation  with  zero  space  charge  was  solved.  However, 
if  the  doping  level  is  of  the  same  order  as  the 
intrinsic  concentration  and  if  significant  electric  and 
magnetic  fields  are  present,  then  the  space  charge 
cannot  be  neglected  and  the  Hall-plate  model  does 
not  apply.  For  example,  Fig.  4  shows  the  computed 
space-charge  distribution  in  a  4  pm  x  4  urn  n-type 


silicon  slab,  where  a  close-to- intrinsic  condition 
is  modeled  by  an  assumed  elevated  temperature  of 
500°K. 


Fig.  3  Equipotential 

lines  and  current  lines 

for  a  square  device; 

L  ■  W  ■  100  u  m,  p-type 

silicon,  doping  level 

1016  cm-3,  B  “10  Tesla, 
z 


Fig.  4  Space  charge 
distribution  for  a 
square  device;  L »  W  ■ 

4  urn,  n-type  silicon, 
doping  level  1033  cm-3, 

B  «  1  Tesla,  T  ■=  500°K. 


In  intrinsic  material,  where  the  net  doping  is 
much  smaller  than  the  intrinsic  concentration, 
both  electron  and  hole  concentrations  are  strongly 
affected  by  the  magnetoconcentration  effect  (see 
Figs.  5b  and  5d) .  This  results  in  an  almost 
exponential  increase  of  the  concentrations  in  x 
direction  along  the  middle  of  the  slab.  On  the 
other  hand,  the  metal  boundaries  impose  a  strong 
concentration  gradient  in  the  y  direction  on  both 
electrons  and  holes.  The  interaction  of  magneto¬ 
concentration  and  boundary  effects  gives  rise  to 
the  unexpected  space-charge  distribution  shown  in 
Fig.  5c.  We  note  that  in  the  intrinsic  case  the 
scaling  rule  breaks  down,  because  the  space 
charge  depends  on  the  local  concentration  gradients. 


In  Fig.  5a,  the  equipotential  lines  show  only 
small  Hall-field  components,  while  the  current 
lines  crowd  on  the  high-concentration  side  of  the 
device.  This  result  suggests  that  a  device  with 
split-metal  contacts,  as  shown  in  Fig.  2,  could 
be  employed  to  measure  a  magnetic  field.  The 
sensitivity  of  such  split-metal-contact  devices  can 
be  defined  by  the  relative  current  imbalance,  viz. 

S  “(Ij  -  I2)  /  (Ij  +  12^"  Table  1  shows  the  sensi¬ 
tivities  computed  for  different  doping  levels  and 
geometries.  The  intrinsic  split-metal-contact 
device  yields  the  highest  sensitivity  of  all  the 
ones  we  have  modeled.  This  intrinsic  device, 
however,  is  rather  sensitive  to  size  reduction; 
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a  4  u m  x  4  urn  split-metal-contact  device  would 
yield  less  than  half  of  the  sensitivity  of  a 
100 u m  x  100  urn  device.  This  is  the  consequence 
of  the  building  up  of  space  charge  mentioned 
earlier  (see  Fig.  5c),  which  acts  against  the 
magnetoconcentration  effect  and  modifies  the 
current  density  and  potential  distributions. 


Fig.  5  Intrinsic  silicon  square  device;  L  -  W  * 

4  pm,  B  “1  Tesla,  T  *  500°K.  a)  Equipotential 
lines  and  current  lines,  b)  Hole  concentration, 
c)  Space-charge  distribution,  d)  Electron  con¬ 
centration. 


Type 

L 

X 

W 

Xl+I2 

Intrinsic 

100 

pm 

X 

100 

pm 

0.443 

Intrinsic 

4 

U  m 

X 

4 

p  m 

0.338 

n 

100 

pm 

X 

100 

p  m 

0.093 

n 

400 

V  m 

X 

100 

pm 

0.101 

n 

100 

pm 

X 

400 

p  m 

0.033 

n  -  1016cm‘3  Vfll  -1.0  Volt 

B  -  1  Tesla  Vfl2  -  1.0  Volt 


Gummel  scheme  to  two  dimensions  and  nonzero  magnetic 
field;  depending  on  the  initial  guess  our  computa¬ 
tions  may  require  only  up  to  one  third  more  CPU¬ 
time  than  the  case  of  zero  magnetic  field. 

In  case  of  Hall  plates,  our  results  agree  with 
prior  results  obtained  by  simple  analytical  models 
or  conformal  mapping  technique.  However,  when 
magnetoconcentration  and  space-charge  effects  occur 
simultaneously,  only  numerical  modeling  will  describe 
the  functioning  of  the  device.  We  have  modeled  devices 
with  split-metal  contacts  which  show  promising  sen¬ 
sitivities;  our  device  principle  has,  however,  to  be 
explored  further  using  more  realistic  ohmic  contacts. 

We  plan  to  apply  our  numerical  algorithm  to 
magnetodiodes,  magnetotransistors ,  and  split-drain  MOS 
magnetic  field  sensors.  Modeling  of  more  general 
directions  of  the  magnetic  field  will,  however,  require 
three  spatial  dimensions.  The  questions  of  optimum 
sensitivity  and  signal-to-noise  ratio  will  be  studied 
for  specific  cases. 
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Table  1  Magnetic  Sensitivity  (1^-12) /(Il+I2^ 
of  silicon  slabs  of  various  dopings  and  geometries 
(see  Fig.  2 ) . 


Discussion  and  Outlook 


A  semiconductor  device  exposed  to  a  magnetic 
field  has  been  modeled  numerically.  The  asymmetries 
introduced  by  the  magnetic  field  into  the  transport 
equations  and  the  boundary  conditions  have  been 
incorporated  successfully  into  a  computer  algorithm 
by  a  generalization  of  the  well-known  Scharfetter- 
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APPENDIX  G 


Matrix  and  Vector  Elements  for  the  Split-Drain  MDSFET 


The  coefficients  of  the  variable  u  associated  with  eqn.  (4.31)  for  the 
bulk  and  eqns.  (4.32)  -  (4.34)  describing  the  floating  boundaries  of 
the  split-drain  MDSFET  (see  Figs.  4.3  -  4.5)  are  described  in  this 
appendix. 

Associated  with  eqn.  (4.30)  are: 


at(m,  N) 


2 


h  (N-l)  [hy(N)  +  hy(N-l)  J 


2 


2 


bt(m,  N) 


h  (M)  h  (M-l) 

X  X 


h  (N)  hy (N-l) 


Ct(M,  N)  = 


2 


hy(N)  [hy(N)  +  hy  (N-l)  ] 


2 


VM'  N)  =  h  (M-l)  111  (M)  +  h  (M-i)] 

x  L  X  x. 
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Eqn.  (4.24)  is  discretised  by  a  central  difference 
scheme,  i.e.. 


At 


u ( 3 ,  N)  -  u ( 1 ,  N)  _  *  _  u ( 2 ,  N+l)  ~  u ( 2 ,  N-l) 

h  (1)  +  h  (2)  yn  z  h  (N)  +  h  (N-l) 

a  x  y  y 


®u (XMAX,  N)  ~  u (XMAX-2 ,  N)  =  *  u(XMAX-l,  N+l)  -  u(XMAX-2) 

;  hx (XMAX-1)  +  hx  (XMAX-2)  yn  z  h  (N)  +  hy(N-l) 


hx(l)  +  hx(2) 
h  (N)  +  h  (N-l) 


Hence,  at 


©• 


Gt(2, 


N)  =  -  w 


n 


B 


Ht(2,  N)  =  -  Gt(2,  N)  , 


It(2,  N)  =  0 


and  Jt(2,  N)  =  1 


At  the  floating  boundary 


G  (XMAX-1,  N) 


hx(XMAX-l)  +  hx (XMAX-2) 
h  (N)  +  h  (N-l) 

y  y 


H  (XMAX-1,  N)  =  -  GT( XMAX-1,  N) 


=  0  and  (XMAX-1,  N)  -  1 


I  (XMAX-1,  N) 


"  \ 

At  the  floating  boundary  ©• 

Using  central  differencing  on  eqn.  (4.23); 


u(M+l,  YMAX-1)  -  u(M-l,  YMAX-1 ) 
h  (M-l)  +  h  (M) 

X  X 


1 


u(Mr  YMAX)  -  u (M,  YMAX- 2) 

h  (YMAX-1)  +  h  (YMAX- 2) 

y  Y 


This  yields  the  coefficients; 


P  (M,  YMAX-1)  = 


-  w 


*  B 
n  z 


hy (YMAX-1)  +  hy(YMAX-2) 
h  (M)  +  h  (M- 1 ) 

A  A 


Qt(M,  YMAX-1)  =  -  Pt(M,  YMAX-1)  , 


R  (M,  YMAX-1)  =  0  and  ST<M,  YMAX-1)  =  1 


APPENDIX  H 


TWO  •  DIMENSIONAL  MODELING  OF  A  MOS  MAGNETIC  FIELD 

SENSOR 


A.  NATHAN.  STUDENT  MEMBER.  IEEE.  L.  ANDOR  . 

H.  P.  BALTES.  MEMBER.  IEEE  ind  H.  G.  SCHMIDT- WEI NMAR 


Abstract  -  We  present  numerical  results  for  the  potential,  current  and 
surface  charge  distributions  as  well  as  the  sensitivity  of  a 
magnetic-field-sensitive  split-drain  N-channel  MOSFET  operating  in  the 
linear  region.  We  also  present  a  suitable  circuit  configuration  that 
incorporates  this  device  and  calculate  the  overall  sensitivity  of  the  resulting 
magnetic-field -sensitive  circuit.  Optimisation  of  the  device  geometry  is 
discussed. 


II.  THEORY 

We  analyse  the  channel  region  of  an  N-channel  MOSFET  in  the 
linear  region,  subject  to  a  magnetic  field  perpendicular  to  the  device  plane. 

B  -  (0.0.BZ  ).  Following  [2],  the  charge  in  the  inversion  layer  is  given 
by 

Qn  (x,y)  -  Cqx  [  V(x,y)  +  2*b  -  VQ  ] 

(1) 

+  '2cqNA  (  V(x.y)  +  2*  ] 


I.  INTRODUCTION 


A  number  of  MOS  magnetic  field  sensors(MFS)  based  on  the 
splitdrain  configuration  have  been  devised  and  realised  recently  (see  [1] 
and  references  therein).  To  the  best  of  our  knowledge,  no  numerical 
modeling  of  such  sensors  has  been  reported  hitherto. 

In  this  note  we  present  the  first  results  of  a  two-dimensional 
numerical  analysis  of  split -drain  NMOS  MFS  operating  in  the  linear 
region,  the  magnetic  field  vector  being  perpendicular  to  the  device  plane. 
Our  numerical  method  allows  to  calculate  the  sensitivity  of  these  MFS  for 
a  variety  of  device  geometries,  i.e.  different  aspect  ratios  and  drain 
separations. 

In  contrast  to  the  MFS  described  previously  [1],  our  model  applies 
to  the  linear  region  of  MOSFET  operation.  Therefore,  in  this  note,  we 
also  present  a  circuit  configuration  suitably  matched  to  the  split -drain 
NMOS  device. 

H  P  Ballcs.  A  Nathan,  and  II  G.  Schmidt-Weinmur  arc  with  the  De¬ 
partment  of  Electrical  Engineering.  University  of  Alberta.  Edmonton. 
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in  the  usual  notation  [2],  While  the  continuity  equations  and  Poisson's 
equation  remain  the  the  same  as  for  zero  field  [3],  the  current  density  has 
to  be  modified  [4],  to  account  for  the  Lorentz  force,  viz., 

1 


J  = 


n  ( (qu  nE  +  qD  ^n) 

1  +  Gj  *  f)  2  n  n 


+  U*  ?  x  (q-^nE  +  qD^n)]  , 


(2) 


in  the  usual  notation  [2].  The  hole  current  and  net  recombination  are 
assumed  zero  in  the  channel.  A  two-dimensional  model  is  obtained  by 
integrating  the  current  density  over  the  channel  depth,  viz.. 


Jnx  <*’*> 


1  +  (U*  B  ) 
n  z 


~2  [  (~  U  ♦  D  f-  ) 

2  n  3x  n  3x 


(3) 


3V  a 

*  “S  Bz  (‘“n  3?  +  Dn  37  11  Qn  (x-y)  * 


Jny  (x'y) 


1  +  (u*  B  )‘ 

n  2 


K-u  f  +D  f  ) 
n  3y  n  3y 


(4) 


~  u*  B  (-u  +  d  —  )  1  n  (x  v) 

n  z  n  3x  n  3x  vn  1  ,y;  ' 
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where  V(x,y)  denotes  the  reverse  bias  voltage  between  a  point  (x,y)  in  the 
channel  and  the  source  electrode  (which  is  grounded).  We  use  a  field 
dependent  drift  mobility  model  proposed  by  [5].  Eqns.  (3)  and  (4)  are 
solved  numerically  using  a  finite  difference  discretisation  scheme  described 
elsewhere  [6].  Since  the  device  is  operated  in  the  linear  region,  we  average 
the  gate  electric  field  over  the  channel  length  and  assume  a  constant 
drain -source  field. 

We  study  the  sensitivity  of  the  split  drain  N -channel  MOS  MFS  in 
terms  of  the  relative  current  imbalance  in  the  two  drains,  (see  Fig.  1)  viz.. 

5  ‘  (ID1  -  V/(1D1  +  V  <5> 

III.  RESULTS  and  DISCUSSION 


Fig.  2.  Relative  current  imbalance.  6  as  a  function  of  magnetic  field.  B 
for  a  square  channel  NMOS,  W  =  L  =  25  pm.  for  different  relative  drain 
separations  (refer  to  text  for  curves  a.  b.  and  c) 


The  current  imbalance.  6  of  the  split-drain  MOSFET  was 

computed  for  various  aspect  ratios  W/L  and  drain  separations,  d  (see  Fig. 

1).  The  MOSFET  is  operated  in  the  linear  region  with 

V  =  V  =  1  V.  V  =  5  V  and  V  =  0  V.  The  oxide  thickness 
D1S  D2S  Co  b 

15  -3 

t  =  100  nm.  bulk  doping  density.  N.  =  10  cm  and  threshold 

ox  A 

voltage.  VT  =  1  V. 


Fig.  1.  MOSFET  channel  geometry 


The  case  of  the  square  device  geometry  with  W  =  L  =  25  pm  is 
considered  in  Fig.  2  for  the  relative  drain  separations 
d/W  =  0.05  (curve  a).  d/W  =  0.2  (curve  b)  and  d/W  =  0.6  (curve  c). 
The  current  imbalance.  <5  is  a  linear  function  of  the  flux  density.  B 


between  -1  and  +1  tesla,  i.e.  the  sensitivity  6  /B  is  constant.  The 
linearity  is  not  affected  by  the  variation  of  d/W.  For  illustration,  the 
distribution  of  potential,  current  and  surface  charge  for  the  case 
d/W  =  0.2  under  1  tesla  field  strength  is  shown  in  Fig.  3. 


Fig.  3.  (a)  Equipotential  and  current  lines  (b)  surface  charge  density 
of  the  NMOS  MFS.  W  =  L  =  25  pm.  d  =  5  pm  operated  at 
VD1S=VD2S  =  1V-VCS  =  5  V  and  Vs  =0V 


The  relative  current  imbalance  6  for  devices  with  aspect  ratios 
W/L  =  5  pm/25  pm  and  W/L  =  100  pm/25  pm  are  shown  in  Figs.  4  and 
5.  The  three  relative  drain  separations  are  same  as  given  in  the  previous 
paragraph. 

We  learn  that  nanowing  the  channel  down  from  25  pm  to  5  pm 
(Fig.  4)  leads  to  only  a  slight  increase  in  sensitivity.  Making  the  channel 
much  broader  (Fig.  5)  results  in  a  substantial  sensitivity  decrease. 

For  the  W/L  >>  1  geometry  (Fig  5).  we  find  a  slight  deviation 


' 


; 
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Fig.  4.  Same  as  Fig.  2  but  W/L  =  5  pm/25  pm 


from  linear  behaviour.  This  could  possibly  explain  some  of  the 
non-linearity  observed  previously  [1]  for  positive  fields  B  0.7  tesla. 


Fig.  5.  Same  as  Fig.  2  but  W/L  =  100  pm/25  pm 


From  the  point  of  view  of  keeping  the  device  small,  but  allowing 
for  the  minimal  drain  separation  required  by  the  design  rules  the  square 
geometry  seems  to  be  favourable. 

IV.  CIRCUIT  DESCRIPTION 

* 

A  circuit  configuration  matched  to  the  square  device 
(W  =  L  =  25  pm,  d  =  5  pm)  is  shown  in  Fig.  6.  The  action  of  the 
Lorenti  force  leads  to  an  increase  of  the  current  in  one  drain  at  the 
expense  of  the  current  in  the  other  drain.  The  currents.  1^  .  ^2  *re 


converted  into  voltages,  VQi  ,  VQ,  by  incorporating  standard  current 
to  voltage  convertors.  The  current  to  voltage  convertors  can  be  realised 
with  inverting  MOS  operational  amplifiers  [7]  using  voltage  shunt 
feedback.  With  reasonably  high  input  impedance.  Ri  and  high  gain.  Ay 
of  the  operational  amplifiers,  we  have  the  amplifier  inputs  effectively 
shorted  i.e.  V  0.  Hence 


'  -  V,  =  -  R  I  , 

01  1  F  D1 

(6) 

r  -  V  =  -  R  I 

02  2  F  D2  * 

(7) 

Here,  Rf  denotes  the  feedback  resistor  (see  Fig.  6).  Moreover,  any 
offsets  at  the  amplifier  inputs,  ID1  £  ID2  with  rero  field,  can  be 
controlled  by  and  V2  .  Rp  is  suitably  chosen  to  match  the  input 
currents  I  D1  , 1 D2  and  to  meet  the  S/N  ratio  requirements. 

Inserting  (6)  and  (7)  into  (5),  we  finally  obtain  the  overall 
sensitivity  of  the  field  sensitive  circuit. 


(V0l  "  V02)  _  {V1  ‘  V 

6/B - - 02 - 1 - 2_  /B 

(V01  +  V02>  -  (V1  +  V 


(8) 


The  predicted  sensitivity  of  the  circuit  considered  above  is 
6  /B  =  0.12  /T.  The  drain  currents  IDl  .  ID2  are  57  pA  and  45  pA 
respectively. 

To  achieve  the  maximum  attainable  sensitivities,  the  drain  voltages 
are  put  at  the  same  potential  in  order  to  avoid  any  shunting  effects 
between  the  drains  and  the  device  is  operated  at  low  gate  voltages,  e.g. 
VG<  5  V. 


- 


no 


V.  CONCLUSION 

In  this  note,  we  reported  results  of  the  first  numerical  modeling  of 
magnetic -field -sensitive  NMOS  circuits.  Our  model  allows  to  calculate  the 
potential,  current  and  surface  charge  distributions  and  the  sensitivity  for  a 
variety  of  MOSFET  geometries.  In  the  case  of  P -channel  MOSFET's,  the 
mobilities  and  in  eqns.  (3)  and  (4)  have  to  be  replaced  by  ^ip 
and  )i*  .  respectively,  etc..  Hence  the  sensitivity  is  reduced  by  a  factor  of 
pn/pp  for  an  equivalent  P -channel  device  in  circuit.  The  device 
characteristics  resulting  from  the  above  numerical  analysis  can  be  readily 
incorporated  into  a  circuit  simulation  program  such  as  SPICE  for 
magnetic -field -sensitive  circuit  simulation. 
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APPENDIX  I 


Listing  of  Split-Drain  MDSFET  Simulation  Program 


IMPLICIT  REAL*8( A-H.L ,0-Z) 

RE AL *8  MUNH , MUN , MUNO , NA , N I 

REAL*8  S ( 50) , SM( 50) , SP ( 50) , SS ( 50) , SSM( 50 ) , 

1  SSP ( 50) . R ( 50) , SR ( 50) ,S0(50) ,SS0(50) , 

2  DSP( 50) , DSM( 50) , DSO( 50) , DS ( 50) , DR ( 50) , 

3  «JN(50.50)  ,QN(  50,50)  ,U(  50,50)  ,  V(  50, 50) 

DIMENSION  AT ( 50 ,50) , BT ( 50 , 50) , CT ( 50 , 50) , DT ( 50 , 50) , 

1  ET ( 50 , 50) , FT ( 50 , 50) , BT I ( 50 , 50) , F ( 50) , DB I ( 50 ,50) 

COMMON/KICSI/  U0( 50 , 50) , V0( 50 , 50) , H , NX , NX  1 , 

1  NY , NY  1 , HX ( 50) , HY( 50) , NSD , NSDS , NSDE 

COMMON/MATR/  AT , BT , CT , DT , ET , FT 

COMMON/PREM/  DB I . BT I , S , SS , SM , SSM . SP , SSP , R . SR , SO , SSO . 

1  DSP , DSM , DSO , DS , DR 

C 

C  Q=ELEMENTARY  CHARGE 

C  CK=BOLTZMANN ' S  CONSTANT 

C  EPS=PERMITIVITY  OF  FREE  SPACE 

C  RR= 1 .92  IS  USED  FOR  HALL  MOBILITY  CALCULATIONS 

C 

0=1 . 6D- 1 9 
CK= 1 . 38D-23 
EPS= 1 2 . D0*8 . 86D- 14 
EPS  1 0  =  3 . 9D0*8 . 86D- 1 4 
RR  = 1 . 92DO 
PRINT34 
C 

C  NA=BULK  DOPING  DENSITY 

C  H=MAGNETIC  INDUCTION  IN  VOLT-SEC  PER  SO.  CM 

C  W  AND  L  ARE  THE  DEVICE  D I  MENS  IONS ( WI DTH  AND  LENGTH) 

C  XWR=RELATIVE  DRAIN  SEPARATI 0N( D/W) 

C  NMC  IS  USED  FOR  TRIPLE-DRAIN  DEVICE  SIMULATIONS 

C  TX  AND  TY  ARE  THE  GRID  SPACING  RATIOS  IN  X  AND  Y 

C  OMEGA=RELAXATI ON  PARAMETER  USED  IN  SLOR 

C  VD,  VS  AND  VG  ARE  DRAIN.  SOURCE  AND  GATE  VOLTAGES 

C  I SLOR=MAX .  NO.  OF  SLOR  ITERATIONS 

C  TOX=OX IDE  THICKNESS 

C  TT  =  TEMPERATURE 

C  INGES=PARAMET£R  FOR  EXTERNAL  INITIAL  GUESS 

C 

34  FORMAT ( '  TYPE  NA , H , W , L . XWR , NMC . TX , TY . OMEGA , VD . VS . VG . ' 
1  /'  ISOR.TOX.TT, INGES ' ) 

READ* , NA . H , W . L . XWR . NMC , TX . TY . OMEGA , VD . VS , VG . 

1  IS. TOX.TT, INGES 

NX  AND  NY  ARE  THE  NO.  OF  GRID  POINTS  IN  X  AND  Y 


NX  =  50 
NY  =  50 
NX  1 =NX -  1 
NY  1 =NY- 1 

CALL  GRID( NX , NY . W . XWR . L . TX . TY . ND 1 E . ND2S . HX . HY ) 

NSDS=ND 1 E+ 1 
NSDE=ND2S- 1 
NSD=ND2S-ND1E 
TETA=Q/(CK*TT) 

COX=OXIDE  CAP.  PER  SO.  CM 

NI  =  I  NT  R I  NS  I C  CARRIER  CONCENTRATION 
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84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 

100 

101 

102 

103 

104 

105 

106 

107 

108 

109 

1  10 

1  1  1 

1  12 

1  13 

1  14 

115 

1  16 

117 

1  18 

1  19 

120 


no  non  non  non  non 
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PSIB=BULK  FERMI-LEVEL 
VT=THRESH0LD  VOLTAGE 

COX=EPSIO/TOX 

NI =3 . 88D 16*TT**  1 . 5*DEXP(-7.D3/TT) 

PSIB= ( 1/TETA)*( DL0G(NA/NI ) ) 

VT=2*PSIB+(2. D0*DSQRT (EPS*Q*NA*PSIB)/COX) 

PRINT  * , VT 

ALPHA=C0X* (2*PSIB-VG) 

BETA=C0X 

GAMMA  =  2*  EPS*NA  *Q 
DELTA=4*EPS*NA*Q*PSIB 

CALL  MBLT V ( VD , VS , L , COX , EPS , VG . PS  I B , NA , MUN . RR . MUNH ) 
PRINT  * , MUN 

TGDN= TANGENT  OF  HALL  ANGLE 
T  GDN  =  MUNH*H 

CALL  INGUES( VD , VS, TETA , ALPHA , BETA , GAMMA , DELTA , INGES) 
CALL  FLBC ( S , SS , SP , SSP , SM, SSM, R , SR , SO, SSO, 

1  DSP.DSM.DSO.DS.DR.TGDN) 

CALL  MATRCO 

CALL  PREMAT ( NX , NX  1 , NY , NY  1 , NSD , NSDS , NSDE ) 

SLOR  ITERATION  BEGINS  FOR  THE  SOLUTION  OF  U 
UAJ=-1 

DO  98  I SOR=  1  ,  IS 
<JA  J  =  -  JA  J 
I C0N=O 

DO  97  NNN  =  2 , NY  1 
N  =  NNN 

IF(JAJ.LE.O)  N=NY 1 +2-NNN 

CALL  FN( UO , NX  1 , NY  1 , NSD , ND 1 E , NSDS . NSDE , ND2S , 

1  N , R , SR , DR , F ) 

CALL  MATCA(NX1 ,N,DBI ,BTI ,ET,U,F) 

DO  96  M  =  2 , NX  1 

IF(DABS(U(M,N)/(UO(M,N) ) -  1 . DO) . GE .  1 . D-5 )  IC0N=1 
UO(M,N)=OMEGA*U(M,N)+( 1 . DO-OMEGA ) *U0( M , N) 

96  CONTINUE 

97  CONTINUE 

BOUNDARY  VALUES  ARE  CALCULATED 

DO  95  N=2 , NY  1 

U0( 1 ,N)=SP(N) *U0( 2 , N+ 1 ) 

1  +SM(N)*U0(2,N-1 )+  S0( N ) *U0( 3 , N) 

95  U0( NX ,N)=SSP(N) *U0( NX  1 , N+ 1 ) 

1  +SSM(N)*U0(NX1 ,N-1 )+  SSO( N ) *U0( NX  1  -  1 , N) 

IF(NSD.EO.O)  GO  TO  98 
DO  94  M=NSDS , NSDE 

94  UO(M. 50)=DSP(M) *U0(M+1 , 49 ) +DSM( M ) *U0( M- 1 ,49) 

1  +DSO(M)*UO(M,48)+DS(M)*UO(M,49)+DR(M) 

98  IF( ICON . EQ . 0)  GO  TO  99 
CALCULATION  OF  THE  POTENTIAL  DIST. 
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c 

99  CALL  NEWRAP(NX. NY, ALPHA. BETA. GAMMA, DELTA. TETA, VO, U0) 

SURFACE  CHARGE  AND  CURRENT  DISTS.  ARE  CALCULATED 

DO  103  N= 1 , NY 
DO  103  M= 1 , NX 

103  QN( M , N ) = - ( ALPHA+BETA* V0( M , N ) +DSQRT ( GAMMA  * VO ( M,N)+DELTA) )/Q 

CALL  CRNT  (  UO  ,  HX  ,  HY  ,  MUN  ,  TGDN  ,  <JN  ,  NX  ,  NY  . 

1  ND 1 E , NSDS .NSDE.ND2S.NSD, NMC , SENS) 


WRITE 

(  10) 

( (ON(M.N) 

,M=1 

.  NX)  , 

N=  1 

.NY) 

WRITE 

(ID 

( (UO(M.N) 

,  M=1 

•  NX)  , 

N=  1 

.NY) 

WRITE 

(12) 

((VO(M.N) 

,  M=  1 

.NX)  , 

N=  1 

.NY) 

WRITE 

(13) 

(  (vJN(M.N) 

.  M=  1 

.NX)  . 

N=  1 

,NY1) 

WRITE 

(  15) 

( HX ( M) ,M= 

1  .NX 

1) 

WRITE 

(15) 

( HY ( N) , N= 

1  ,  NY 

D 

WRITE 

(20) 

ND1E.ND2S 

.MUN 

.TGDN 

WRITE 

(20) 

NX  .NY 

PRINT* . SENS 

PRINT  * , ND 1 E , NSD , ND2S 

STOP 

END 

GRID  GENERATION 

SUBROUTINE  GR I D ( NX . NY , W , XWR . L , TX , T Y , NW 1 , NW2 , HX , HY ) 
IMPLICIT  REAL*8(A-H,L,0-Z) 

REAL*8  HY ( 50 ) , HX ( 50 ) ,L,W1 ,W2,W3 

W2=XWR*W 
W 1 = { W-W2 ) /2 . DO 
W3  =  W  1 
NX  1 =NX- 1 
NY  1 =NY -  1 

CONST  =  NX/ (DSQRT(W1 ) +DSQRT ( W2 ) +DSQRT ( W3 ) ) 

WW 1 = CONST *D SORT ( W 1 ) 

WW2  =  CONST*DSORT ( W2 ) 

WW3=C0NST*DSQRT(W3) 

NW 1 = WW 1 +0 . 7 
NW2=WW2+WW 1+0.7 
NW3=WW 1 +WW2+WW3+0 . 7 

NS=  1 

CALL  NUGR I D ( NS , NW 1 , W 1 , TX , HX ) 

NS=NW 1 

CALL  NUGR ID ( NS , NW2 , W2 , TX , HX ) 

NS=NW2 

CALL  NUGR ID ( NS , NW3 , W3 . TX , HX ) 

NS=  1 

CALL  NUGR ID(NS.NY,L.TY,HY) 

RETURN 

END 
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NON-UNIFORM  GRID 

SUBROUTINE  NUGR  ID ( NS , N, XY , TXY , HXY ) 

IMPLICIT  REAL*8(A-H,L,0-Z) 

REAL *8  HXY ( 50) , XY 
N 1  =N-  1 
N2=N- ( NS+ 1 ) 

SC A =4 . DO* ( TXY- 1 )/(N2)**2 
HXYS=0 . DO 
DO  4  M=NS , N 1 

HXY ( M ) = 1 . DO+SCA* (N1-M)*(M-NS) 

4  HXYS  =  HX YS  +  HX Y ( M ) 

SCA=XY/HXYS 
DO  5  M=NS , N 1 

5  HX Y ( M ) =SCA  *HXY ( M) 

RETURN 

END 

CALC.  OF  THE  FIELD  DEPENDENT  MOBILITY 
MUN  =  DR I  FT  MOBILITY  FOR  ELECTRONS 
MUNH=HALL  MOBILITY  FOR  ELECTRONS 

SUBROUTINE  MBLT Y ( VD , VS , L , COX , EPS , VG , PS  I B , NA , MUN . RR , MUNH ) 
IMPLICIT  REAL*8( A-H.L ,0-Z) 

REAL  *8  MUNO.NA, MUN, MUNH 

MUNO= 1 . 6D3 

EP=(VD-VS)/L 

EPA=EP/3 . 5D3 

EPB=( EP/7 . 4D3 ) *  *2 

ETR=(C0X/EPS)*(VG-PSIB-VD/2) 

FEP=  t . DO/DSORT ( 1 . DO+NA/ ( NA/3 . 5D2  +  3 . D 1 6 )  + 

1  (EPA**2)/(EPA+8. 8D0) +EPB ) 

GETR= 1 . DO/DSORT ( 1 .DO+1 ,539D-5*ETR) 

MUN  =  MUNO*  F  EP*GETR 
MUNH=MUN*RR 

RETURN 

END 

GENERATION  OF  INITIAL  GUESS  FOR  VARIABLES  U  AND  V 

SUBROUTINE  I NGUES ( VD , VS , T ETA , ALPHA , BETA , GAMMA , DELTA , INGES) 
IMPLICIT  REAL*8(A-H,L ,0-Z) 

REAL *8  GUESS( 50, 50) 

COMMON/KICSI/  U0(  50, 50) , V0(50,50) , H.NX.NX1 , 

1  NY , NY  1 , HX ( 50) , HY ( 50) , NSD , NSDS , NSDE 

READ  (14)  ( ( GUESS (M,N) , M= 1 .NX) ,N=1 ,NY) 

DO  1  N= 1 , NY 
DO  1  M= 1 , NX 

UD=ALPHA/TETA  + ( BET A/TE T A - AL PHA ) * VD 

1  -BETA  * ( VD*  *2 ) /2 

2  +DSQRT ( GAMMA* VD+DELTA ) /TETA 

3  -2* ( ( GAMMA* VD+DELT  A ) *  *  1 . 5 ) / ( 3 *GAMMA ) 

US=ALPHA/TETA  +( BETA/TETA- ALPHA )* VS 
1  -BETA*(VS**2)/2 
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2  +DSQRT(GAMMA*VS+DELTA)/TETA 

3  -2*( (GAMMA*VS+DELTA) ** 1 . 5 ) / ( 3*GAMMA ) 

U0(M,N)=( UD-US )/NY1*(N-1)+US 
IF( INGES.NE . 1 )  GO  TO  1 
U0( M , N) =GUESS ( M,N) 

1  V0(M,N)=(VD-VS)/NY1*(N-1 )+VS 

RETURN 

END 

EVALUATION  THE  COEFFS.  FOR  FLOATING  BOUNDARY  MESHPOINTS 

SUBROUTINE  FLBC ( S , SS , SP , SSP , SM , SSM , R , SR , SO , SSO . 

1  DSP , DSM , DSO , DS , DR , T  GDN ) 

IMPLICIT  REAL*8(A-H,L,0-Z) 

COMMON/KICSI/  U0( 50, 50) , V0( 50 , 50 ) , H , NX  ,  NX  1  , 

1  NY , NY  1 , HX ( 50) , HY ( 50) , NSD , NSDS , NSDE 

REAL*8  S( 50) , SM( 50) . SP ( 50) . SS ( 50) ,SSM(50) .SSP (50) ,R(50) , 

1  SR (50) . S0( 50) , SSO( 50) .DSP (50) ,DSM(50) .DS0(50) . 

2  DS ( 50) , DR ( 50) 

DO  1  N=2 , NY  1 

R ( N ) =0 . DO 
S ( N ) =0 . DO 

SM(N)  =  (T  GDN* ( HX ( 1 ) +HX (2) ) )/ (HY(N)+HY(N-1 ) ) 

SP(N)=-SM(N) 

S0( N ) = 1 .DO 

NX2=NX 1-1 
SR ( N ) =0 . DO 
SS ( N ) =0 . DO 

SSM(N)=-(T GDN* ( HX ( NX  1 )+HX(NX2) ) ) / ( HY ( N ) +HY ( N- 1 ) ) 
SSP(N)=-SSM(N) 

1  SSO( N) = 1 . DO 

IF(NSD.EQ.O)  GO  TO  3 
NY2  =  NY 1  -  1 

DO  2  M=NSDS , NSDE 
DR ( M ) =0 . DO 
DS ( M ) =0 . DO 

DSM(M)=TGDN*(HY(NY1 )+HY(NY2) ) / ( HX ( M ) +HX ( M- 1 ) ) 
DSP(M)=-DSM(M) 

2  DSO( M ) = 1 . DO 

3  RETURN 
END 

COMPUTING  OF  THE  COEFFS.  IN  THE  DEVICE  INTERIOR 

SUBROUTINE  MATRCO 
IMPLICIT  REAL*8(A-H.L ,0-Z) 

DIMENSION  AT (50. 50) , BT( 50,50 ) ,CT( 50,50 ) . 

1  DT(50, 50) , ET( 50, 50) , FT( 50, 50) 

COMMON/KICSI/  U0( 50 , 50) ,V0(50,50) , H , NX , NX  1 , 

1  NY , NY  1 , HX ( 50) , HY ( 50 ) , NSD , NSDS , NSDE 

COMMON/MAT R/  AT , BT , CT , DT , ET . FT 
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DO  1  N  =  2 , NY  1 
DO  1  M=2 , NX  1 

AT(M,N)=2 . DO/HY ( N  -  1 )/(HY(N)+HY(N- 1 ) ) 
CT(M,N)=AT(M,N)*HY(N-1)/HY(N) 

DT(M,N)=2 . DO/HX ( M- 1)/(HX(M) +HX ( M  -  1 ) ) 

ET(M,N)=DT(M,N) *HX (M-1 )/HX(M) 
BT(M,N)=-(AT(M,N)+CT(M,N)+DT(M,N)+ET(M,N) ) 

1  FT ( M , N) =0 . DO 

RETURN 

END 

INVERSION  OF  THE  COEFFS.  B  PRIOR  TO  SLOR(SEE  [5]) 

SUBROUTINE  PREMAT ( NX , NX  1 , NY , NY  1 , NSD , NSDS , NSDE ) 
IMPLICIT  REAL*8(A-H,L ,0-Z) 

REAL *8  S ( 50) , SM( 50) , SP( 50) , SS( 50) ,SSM(50) , 

1  SSP ( 50) , R( 50) , SR( 50) , SO( 50) ,SS0(50) , 

2  DSP (50) , DSM( 50 ) ,DS0(50) , DS ( 50) . DR ( 50) 
DIMENSION  AT (50, 50) ,BT( 50.50 ) ,CT( 50,50 ) . 

1  DT(50, 50) , ET(50, 50) , FT(50, 50) , 

2  BTI(50,50) ,DBI(50, 50) 

COMMON/M ATR/  AT , BT , CT , DT , ET , FT 

COMMON/PREM/  DB I , BT I , S , SS . SM , SSM , SP , SSP , R . SR . SO , SSO , 
1  DSP,DSM,DSO,DS,DR 


DO  1  N  =  2 , NY  1 

CE=DT(2,N)*S0(N) 

CA=DT(2,N)*SM(N) 

CB=DT(2,N)*S(N) 

CC=DT(2,N)*SP(N) 

ET(2,N)=ET(2,N)+CE 

AT(2.N)=AT(2,N)+CA 

BT(2.N)=BT(2,N)+CB 

CT(2,N)=CT(2,N)+CC 

CD  =  ET ( NX  1 ,N) *SSO(N) 

CA  =  ET ( NX  1 ,N)*SSM(N) 
CB=ET(NX1 ,N)*SS(N) 

CC  =  ET ( NX  1 , N) *SSP(N) 

DT ( NX  1 , N ) =DT ( NX  1 ,N)+CD 
AT ( NX  1 ,N)=AT(NX1 , N ) +CA 
BT ( NX  1 , N ) =BT ( NX  1 , N) +CB 
CT ( NX  1 , N) =CT ( NX  1 , N ) +CC 

1  CONTINUE 

IF(NSD.EO.O)  GO  TO  3 

DO  2  M=NSDS , NSDE 

DA  =  CT ( M , NY  1 )*DSO(M) 
DB=CT(M,NY 1 )*DS(M) 
DD=CT(M,NY 1 )*DSM(M) 
DE=CT(M,NY 1 )*DSP(M) 


AT(M,NY 1 )=AT(M,NY 1 )+DA 
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BT(M,NY1 )=BT(M,NY 1 )+DB 
DT(M,NY1)=DT(M,NY 1 )+DD 
ET(M,NY1)=ET(M,NY1)+DE 

2  CONTINUE 

3  CONTINUE 

DO  6  N  =  2 , NY  1 

BTI ( 2 , N) = 1 . DO/BT ( 2 , N ) 

DO  6  M  =  3 , NX  1 

DBI(M,N)=DT(M,N)*BTI(M-1  ,N) 
BT(M,N)=BT(M,N)-(DBI(M,N)*ET(M  -  1  ,N)  ) 

BT I ( M , N ) = 1 . DO/BT ( M , N ) 

6  CONTINUE 

RETURN 

END 

C 

C  CALCULATION  OF  RHS  OF  SLOR  EON. 

C 

SUBROUTINE  FN(UO, NX  1 , NY  1 , NSD , ND 1 E , NSDS , NSDE , ND2S  , 
1  N.R.SR.DR.F) 

REAL *8  U0( 50 ,50) , R ( 50) , SR ( 50) , DR ( 50 ) , 

AT ( 50 , 50) , CT ( 50 . 50) , BT ( 50 , 50 ) . 

DT ( 50 , 50) , ET ( 50 . 50 ) , FT ( 50 , 50 ) . 

F ( 50 ) , FD,FE,FA,FC,DFC 
COMMON/MATR/  AT , BT . CT , DT , ET , FT 

IF(NSD.NE.O)  GO  TO  4 
DO  1  M  =  2 , NX  1 

FA=AT(M,N)*UO(M,N- 1 ) 

FC=CT(M,N) *U0( M , N+ 1 ) 

F(M)=FT(M,N)-FA-FC 
1  CONTINUE 

FD=DT(2,N)*R(N) 

FE=ET (NX  1 , N ) *SR ( N ) 

F(2)=F(2)-FD 
F ( NX  1  )  =  F  ( NX  1 ) - F  E 

RETURN 

4  DO  6  M=2 , ND 1 E 

FA=AT(M,N) *U0( M , N- 1 ) 

FC=CT(M.N) *U0( M . N+ 1 ) 

F(M)=FT(M,N)-FA-FC 

6  CONTINUE 

DO  5  M  =  ND2S , NX  1 

FA=AT(M,N)*UO(M,N- 1 ) 

FC=CT(M,N)*U0(M,N+1 ) 

F(M)=FT(M.N)-FA-FC 

5  CONTINUE 
FD=DT(2,N)*R(N) 
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FE=ET(NX1 , N ) *SR ( N) 

F ( 2 ) =F ( 2 ) -FD 
F (NX  1 ) =F(NX 1 ) -FE 

DO  2  M=NSDS , NSDE 
IF(N.NE.NYI)  GO  TO  7 

DFC=CT(M,N)*DR(M) 

FA=AT(M,N) *U0( M , N- 1 ) 

F(M)=FT(M,N)-FA-DFC 
GO  TO  2 

7  FA=AT(M.N) *U0( M , N- 1 ) 

FC=CT(M.N)*U0(M.N+1) 

F(M)=FT(M.N)-FA-FC 

2  CONTINUE 
RETURN 
END 

SOLUTION  OF  U  AT  EACH  MESHPOINT 

SUBROUTINE  MATCA (NX  1 , N , DB I , BT I , ET , U , F ) 

REAL *8  AF ,CW,DBI ( 50 , 50) , BT I ( 50 , 50 ) , E T ( 50 , 50) , 

1  U( 50, 50) , F ( 50) 

DO  1  M=3 , NX  1 

AF=DBI (M,N) *F(M- 1 ) 

F(M)=F(M)-AF 

1  CONTINUE 

U ( NX  1  , N ) =BT I ( NX  1 , N ) *  F ( NX  1 ) 

DO  6  K  =  3 , NX  1 
M=NX 1 -K+2 

CW=ET(M,N)*U(M+1 , N ) 

U(M,N)=BTI (M,N)*(F(M)-CW) 

6  CONTINUE 
RETURN 
END 

CALC.  OF  V  USING  NE WTON-RAPHSON 

SUBROUTINE  NEWRAP ( NX . NY . ALPHA , BET A , GAMMA . DELTA , TETA , VO . UO) 
IMPLICIT  REAL*8(A-H,L,0~Z) 

REAL*8  V0(50. 50) ,U0(50, 50) 

DO  1  N= 1 , NY 
DO  1  M= 1 , NX 
DO  2  INR= 1 , 100 

I F ( V0( M , N) . L  E . 0 . DO)  V0(M.N)=0.D0 

FV=UO(M,N)-ALPHA/TETA  - ( BETA/TETA -ALPHA ) *V0( M . N ) 

1  +BETA*(VO(M,N)**2)/2 

2  -DSQRT ( GAMMA  * V0( M , N ) +DELT A) /TETA 

3  +2* ( ( GAMMA  * VO ( M,N)+DELTA)** 1 . 5 ) / ( 3*GAMMA ) 
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FPV=ALPHA-BETA/TETA 
1  +BETA*V0(M,N) 

-GAMMA/ ( DSQRT (GAMMA* VO (M,N)+DELT A) )/TETA/2 
+DSQRT ( GAMMA  * VO ( M ,N)+DELTA) 

IF(DABS(FV)  .LE.  1 . D  -  1 9 )  GO  TO  1 
2  VO(M,N)=VO(M,N)-FV/FPV 

I F ( V0( M , N) . LE . 0 . DO)  VO(M,N)=O.DO 

1  CONTINUE 
RETURN 
END 
C 

C  CALC.  OF  THE  CURRENT  DI5T.  IN  THE  CHANNEL, 

C  THE  DRAIN  CURRENTS  AND  EVALUATION  OF  THE 

C  DEVICE  SENSITIVITY  (DEFINED  IN  TERMS  OF 

C  THE  CURRENT  IMBALANCE  IN  THE  DRAINS). 

C  A  CHECK  FOR  THE  CONSERVATION  OF  TERMINAL 

C  CURRENTS  IS  MADE. 

C 

SUBROUT  I NE  CRNT ( UO . HX , HY , MUN . TGDN , UN . NX , NY , 

1  ND 1 E , NSDS , NSDE , ND2S , NSD , NMC , SENS ) 

IMPLICIT  REAL*8( A-H, L , 0-Z) 

REAL *8  U0( 50 , 50) , HX ( 50) , HY ( 50 ) . MUN , UN( 50 , 50) 

NX  1 =NX- 1 
NY  1 =NY -  1 

MS  =  2 
ML=NX 

DO  1  N=1  , NY  1 

1  CALL  UNY ( UN , UO , HX ,HY,N,MS.ML, MUN , T  GDN ) 

DO  2  N= 1 , NY  1 
DO  2  M= 1 .NX 

2  UN(M.N)=UN(M,N)/UN(NX ,N) 


MS  =  2 
ML  =NX 
N=  1 

CALL  YCRNT ( CRNT , UO , HX , HY ,N,MS,ML , MUN, TGDN) 
SCR=CRNT 


N  =  NY  1 
MS  =  2 
ML  =ND 1 E 

CALL  YCRNT (CRNT , UO . HX , HY , N , MS . ML , MUN . T GDN ) 
DCR 1 =CRNT 

MS  =ND2S+ 1 
ML  =NX 

CALL  YCRNT ( CRNT , UO . HX , HY , N , MS , ML , MUN .TGDN ) 
DCR2=CRNT 

CHE CK= SCR -DCR 1 -DCR 2 

PRINT* , SCR. DCR 1 ,DCR2, CHECK 

SENS=(DCR1-DCR2)/(DCR1+DCR2) 

RETURN 

END 
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SUBROUT  I NE  JNY ( JN . UO , HX , HY , N , MS , ML , MUN , TGDN ) 

IMPLICIT  REAL*8(A-H,L,0-Z) 

REAL *8  U0( 50 , 50) , HX ( 50 ) ,HY(50) ,JN( 50 , 50) .MUN 

JN( 1 ,N)=O.DO 
DO  20  M=MS , ML 

AUO= ( U0( M- 1 , N) +U0( M,N) )/2 . DO 
AU 1 = ( UO( M- 1 ,N+1 )+U0( M , N+ 1 ) )/2 . DO 
DYU=(AU1 - AUO ) /HY ( N ) 

AUO= ( U0( M- 1 , N ) +  U0( M- 1 , N+ 1 ) )/2 .DO 
AU1=(U0(M,N)+U0(M,N+1 ) )/2 .DO 
DXU=(AU1 - AUO) /HX ( M- 1 ) 

20  JN(M,N)=UN(M- 1 , N ) + ( DYU+TGDN*DXU ) *MUN*HX ( M- 1 ) / ( 1 ,D0+TGDN**2) 

RETURN 

END 


SUBROUT  I NE  YCRNT ( CRNT , UO . HX . HY . N . MS . ML , MUN . TGDN ) 
IMPLICIT  REAL*8(A-H,L,0-Z) 

REAL *8  U0( 50 , 50 ) , HX ( 50 ) , HY ( 50) .MUN 

CRNT  =0 . DO 
DO  20  M=MS . ML 

AUO= ( U0( M- 1 ,N)+UO(M,N) )/2 .DO 
AU 1 = ( U0( M- 1 ,N+1 )+U0(M,N+1 ) )/2 .DO 
DYU= ( AU 1 -AUO) /HY(N) 

AUO= ( U0( M- 1 , N ) +U0( M- 1 ,N+ 1 ) )/2 . DO 
AU1=(U0(M,N)+U0(M.N+1 ) )/2.D0 
DXU= ( AU 1  - AUO) /HX ( M- 1 ) 

20  CRNT  =  CRNT  + ( DYU+TGDN*DXU ) *MUNJ  HX ( M- 1 )/( 1  ,D0+TGDN^*2) 

RETURN 

END 


